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Introduction 


This  project  is  aimed  at  exploring  energy-  and  intensity-modulated  electron  beams  for  breast  cancer 
treatment  to  deliver  optimized  conformal  radiotherapy  dose  distributions  that  closely  match  the  target 
volume  and  minimize  the  dose  to  critical  normal  structures.  We  have  proposed  to  work  on  the 
following  tasks:  (1)  to  characterize  electron  beams  from  Helium-filled  accelerators,  (2)  to  develop 
optimization  algorithms  for  energy-  and  intensity-modulated  radiotherapy  (EIMRT)  using  these 
electron  beams,  (3)  to  verify  these  optimized  dose  distributions  using  the  Monte  Carlo  simulation 
technique,  and  (4)  to  compare  the  optimized  dose  plans  obtained  by  EIMRT  with  conventional 
treatment  plans  and  those  obtained  by  photon  intensity-modulated  radiotherapy  (IMRT). 

Body 

Although  photon  beams  have  been  an  effective  modality  for  breast  cancer  treatment  in  radiation 
therapy  the  following  problems  (or  potential  areas  of  improvement)  remain:  (1)  the  inclusion  of  the 
lung  and  sometimes  of  a  small  volume  of  the  heart  in  the  high-dose  volume  due  to  tumor  location, 
patient  size  or  in  the  case  of  chest-wall  treatments;  (2)  lower  dose  near  the  skin  surface  due  to  lack  of 
electron  build-up  in  a  photon  beam;  and  (3)  high  exit  or  scatter  dose  to  the  normal  structures  such  as 
the  lung  and  heart,  and  more  importantly  the  contralateral  breast,  which  may  be  a  major  cause  for  the 
occurrence  of  secondary  cancer  in  the  contralateral  breast  for  women  under  the  age  of  45. 

Recent  development  of  computer-controlled  medical  linear  accelerators  along  with  improved  treatment 
planning  techniques,  may  provide  new  solutions  in  the  delivery  and  control  of  external  beam  radiation 
through  beam-intensity  modulation''^.  It  is  expected  that  using  photon  IMRT,  the  problem  (1)  above 
may  be  significantly  improved  but  (2)  will  remain  and  (3)  may  become  more  serious  as  treatment  time 
increases  with  the  number  of  fields  used  (increased  leakage  or  scattering  dose).  With  modulated 
electron  beams'®''^,  on  the  other  hand,  problem  (1)  may  also  be  significantly  improved  and  problems 
(2)  and  (3)  can  be  completely  eliminated  due  to  the  nature  of  electron  beams. 

In  this  project  we  proposed  to  study  He-filled  accelerators  equipped  with  computer-controlled  multi¬ 
leaf  collimators  and  the  Monte  Carlo  treatment  planning  technique  for  energy-  and  intensity-modulated 
electron  beams  for  breast  cancer  to  significantly  improve  the  dose  uniformity  in  the  target  volume,  to 
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exclude  the  lung  and  heart  from  the  high  or  moderate  dose  volume  and  to  eliminate  the  scatter  dose  to 
the  lung,  the  heart  and  the  contralateral  breast  to  reduce  complications  and  late  effects  associated  with 
breast  cancer  radiotherapy.  The  purpose  of  this  study  is  to  determine  the  degree  of  feasibility  of  using 
electron  EIMRT  to  improve  breast  cancer  treatment. 

This  project  has  4  specific  aims:  (1)  To  perform  Monte  Carlo  beam  simulation;  (2)  To  develop  an 
optimization  algorithm  for  EIMCRT;  (3)  To  perform  Monte  Carlo  dose  verifications;  and  (4)  To 
evaluate  the  optimized  treatment  plans. 

We  report  on  the  research  accomplishments  associated  with  the  tasks  outlined  in  the  approved 
“Statement  of  Work”  below  for  our  research  between  Oct.  1, 1999  and  Sept.  30, 2000: 

1.  Install  optimization  software  for  clinical  treatment  planning  system 

We  have  developed  a  software  interface  called  MCCALC  to  perform  treatment  planning  optimization 
for  clinical  applications.  The  interface  provides  a  graphics  user  interface  (GUI)  for  the  user  to  input 
optimization  parameters  and  to  acquire  patient  CT  and  contour  information  from  a  clinical  treatment 
planning  system  (FOCUS,  CMS,  Inc).  The  computation  of  the  beamlet  dose  distributions  will  be 
performed  on  a  networked  parallel  computer  system  using  the  Monte  Carlo  method  that  we  have 
established  based  on  previous  work’^'^^  and  the  results  will  be  used  for  the  optimization  process.  The 
final  intensity  maps  for  each  treatment  field  will  be  sent  back  to  the  clinical  treatment  planning  system 
to  beam  delivery  using  custom  made  electron  cutouts  or  using  computer  controlled  electron  multileaf 
collimator  (eMLC).  More  detailed  descriptions  of  the  optimization  algorithm  and  its  implementation 
have  been  reported  in  a  publication^'*  and  here  we  only  describe  briefly  for  completion. 

A  fast  iterative  optimization  algorithm  has  been  developed  for  electron  beam  optimization.  The 
algorithm  has  been  tested  for  both  electron  treatment  planning  optimization  and  photon  beam 
optimization.  The  beamlet  dimensions  will  be  specified  through  the  RTF  module,  based  on  the  MLC 
leaf  width.  For  the  Varian  MLC,  the  beamlet  size  will  be  variable  between  10mm  x  1mm  and  10mm  x 
10mm.  The  beamlet  profile  will  be  different  depending  on  the  location  of  the  beamlet  (its  slight  spatial 
dependence  on  the  accelerator  head  scatter,  beam  flatness  and  SSD  variation).  This  will  be 
automatically  included  when  using  the  simulated  phase-space  data  as  source  input.  The  MCDOSE 
code  (see  below)  has  been  modified  for  the  dose  calculations  in  the  patient.  The  beam  incident 
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angles  and  patient  CT  data  will  be  obtained  from  the  RTP  system.  Each  beamlet  will  require  1  -  10 
million  phase  space  photons  depending  on  the  dimensions  of  the  beamlet.  For  a  typical  IMRT  case 
with  9  gantry  angles,  the  beamlet  calculation  (for  1-10  thousand  beamlets)  will  be  completed  in  a  few 
minutes  on  the  PC  network  after  the  variance  reduction  techniques  are  implemented. 

The  optimization  procedure  consists  of  the  following  two  stages.  First,  the  planner  inputs  the  patient 
geometry  and  defines  the  treatment  setup,  such  as  the  beam  energy,  number  and  orientations  of  beams, 
etc.  The  target  volume  and  the  critical  structures  are  defined  by  the  clinician.  The  planner  also 
determines  the  size  of  the  beamlets  and  number  of  dose  or  constraint  points  placed  inside  the  target  or 
critical  structures,  according  to  the  patient  anatomy,  the  required  computation  accuracy  and  the 
available  computation  time.  Each  broad  beam  is  divided  into  beamlets  and  the  dose  and  constraint 
points  are  uniformly  and  randomly  placed  inside  the  corresponding  area.  Then  a  reference  monitor  unit 
is  assigned  to  each  open  rectangular  beam  and  the  dose  deposition  coefficients,  which  is  defined  as  the 
dose  contribution  from  a  beamlet  to  a  point,  are  calculated  using  MCDOSE  and  a  conventional  dose 
calculation  model  for  comparison.  Second,  using  the  calculated  dose  deposition  coefficients  as  input, 
the  optimal  intensity  profile  for  each  beam  is  achieved  by  the  optimization  process  using  an  objective 
function.  The  results  of  the  optimization  process  are  the  intensity  profiles  for  the  individual  gantry 
angles  (photon  fields). 

The  results  of  the  optimization  calculation  are  the  beam  intensity  maps  at  each  gantry  angle  and  then 
they  can  be  used  to  generate  leaf-setting  sequences  using  a  leaf  sequencing  algorithm.  The  final  dose 
calculation  will  be  performed  again  using  MCDOSE.  The  effect  associated  with  leaf  and  jaw 
movement  will  be  accounted  for  in  these  calculations.  We  have  been  evaluating  different  leaf 
sequencing  algorithms  suitable  for  both  “stop  and  shoot”  and  dynamic  delivery.  A  new  algorithm  also 
synchronizes  the  leaf  sequences  to  remove  the  “tongue  and  groove”  effect.  The  results  showed  that  the 
difference  in  the  beam  delivery  time  using  a  dynamic  MLC  between  “stop  and  shoot”  (including  beam- 
off  time  for  leaf  movement)  and  dynamic  delivery  was  clinically  insignificant.  We  will  install  a  leaf 
sequence  algorithm  for  our  project  which  uses  the  stop  and  shoot  algorithm  and  also  synchronizes  the 
leafs  movement  to  remove  the  “tongue  and  groove”  effect.  We  will  further  work  on  other  MLCs  when 
the  leaf  sequence  file  format  becomes  available  for  operation  with  electron  beams. 

2.  Investigate  suitable  objective  functions  for  electron  beam  optimization 
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We  have  installed  several  objective  functions  for  the  optimizer  used  for  electron  beam  EIMRT.  We 
have  selected  a  dose  based  objective  function  with  soft  constraints  using  dose  volume  histogram 
information.  With  this  final  selection  decided,  the  optimization  can  be  performed  by  using  a  dose 
based  penalty  function  method  and  the  center-of-mass  method  to  minimize  the  augmented  objective 
function.  An  improved  conjugate  vector  method  can  also  be  used.  For  the  target  area,  a  quadratic  form 
of  objective  function  is  specified.  In  addition,  two  target  dose-uniformity  constraints  are  used  to  ensure 
a  uniform  target  dose  distribution  and  to  distinguish  the  clinical  importance  of  cold  and  hot  spots.  For 
the  critical  structures,  maximum-dose  constraint  and  several  levels  of  dose-volume  constraints  are 
assigned  to  each  structure.  For  each  objective  function  and  constraint,  an  importance  weight  relative  to 
the  target  objective  function  is  assigned.  All  the  constraints  are  mathematically  transformed  to  the 
penalty  functions  of  quadratic  forms.  The  augmented  objective  function,  which  should  be  minimized, 
is  a  combination  of  the  original  objective  functions  and  all  penalty  functions.  The  detailed  methods 
have  been  reported  in  our  publications^"^  and  the  applications  have  been  presented. 

3.  Verify  electron  EIMRT  using  EGS4/MCDOSE 

We  have  developed  a  Monte  Carlo  EGS4'^  user  code  MCDOSE^^  for  electron  beamlet  and  treatment 
plan  dose  calculations.  Good  agreement  was  achieved  between  the  MCDOSE  results  and 
measurements^"*'^*.  Features  of  MCDOSE  include  a  multiple-source  model  to  reconstruct  the  beam 
phase  space^^,  inclusion  of  beam  modifiers  such  as  jaws,  wedges,  blocks,  compensators  and  electron 
cutouts  in  the  patient  simulation,  the  implementation  of  several  variance  reduction  techniques,  and 
suitable  for  both  conventional  and  intensity  modulated  radiation  therapy  (IMRT)  treatment  planing. 
Before  MCDOSE  is  used  reliably  for  dose  calculation  in  clinic,  it  must  be  properly  validated.  The 
clinical  validations  for  beam  modifiers  and  dose  calculation  are  presented.  A  comparison  of  the  dose 
distribution  with  45-degree  wedge  in  a  water  phantom  made  between  MCDOSE  and 
EGS4/BEAM/DOSXYZ*’’'^’^“  demonstrates  that  MCDOSE  can  give  accurate  result  with  wedge  (Fig. 
1).  The  dose  distributions  for  a  blocked  10x10  15MV  photon  beam  in  a  water  phantom  are  also 
compared  between  them.  They  agree  very  well.  The  comparison  of  Electron  cutout  factors  between 
MCDOSE  and  measurements  show  a  good  agreement  too.  After  applying  variance  reduction 
techniques  in  MCDOSE,  the  agreement  of  dose  distributions  in  specifically  designed  inhomogeneous 
phantoms  between  MCDOSE  and  DOSXYZ  is  within  the  statistical  uncertainty  of  0.5%  (Figs.  2-4). 
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All  these  results  demonstrate  that  MCDOSE  is  accurate  for  routine  dose  calculation  in  radiotherapy 
treatment  planing.  The  heterogeneity  correction  factors  calculated  by  MCDOSE  for  layered-lung  or 
layered-bone  phantoms  were  consistent  with  results  from  measurement  to  within  1%.  Due  to  the 
elegant  variance  reduction  techniques,  MCDOSE  is  also  faster  than  EGS4/DOSXYZ^*^  dose  calculation 
by  a  factor  of  up  to  30.  A  nine  field  IMRT  planing  can  be  done  in  1-4  hours  on  a  personal  computer, 
including  pre-  and  post-optimization  dose  calculation  . 

Initial  schedule  included  the  implementation  of  the  PEREGRINE  code  system  for  this  study.  This  was 
not  done  because  of  the  delay  of  the  commercial  availability  of  the  PEREGRINE  system  and  the 
complication  of  the  commercialization  of  the  PEREGRINE  system.  The  system  has  been  exclusively 
licensed  to  the  NOMOS,  Corp.  and  it  is  no  longer  available  as  a  free  software.  Instead,  we  have 
developed  the  MCDOSE  system,  which  is  based  on  the  EGS4  system,  which  is  free  for  research  and 
education  use.  The  MCDOSE  code  has  been  validated  against  the  DOSXYZ  code  and  proved  to  be 
equally  accurate  but  about  20  times  faster  in  computation  speed.  Therefore  we  have  replaced  the 
PEREGRINE  software  with  MCDOSE  for  this  study. 

4.  Verify  photon  IMRT  using  EGS4/MCDOSE 

To  verify  the  IMRT  dose  distributions  calculated  by  the  commercial  treatment  optimization  system 
(CORVUS,  NOMOS  Corp.)  we  have  modified  the  MCDOSE  code  to  compute  the  dose  from  the  leaf 
sequence  files  from  the  commercial  system.  First  we  performed  calculations  for  a  homogeneous 
PMMA  phantom.  Since  CORVUS  calculates  dose  to  water  while  PMMA  is  not  water  equivalent  it  is 
important  to  establish  a  conversion  scheme  so  that  the  measured  ionization  in  the  PMMA  phantom  can 
be  converted  to  dose  to  water  in  order  to  verify  the  dose  distributions  computed  by  CORVUS.  The 
Monte  Carlo  method  computes  both  dose  to  water  and  dose  to  PMMA  and  the  results  can  be  used  to 
study  the  energy  fluence  perturbation  correction  factor  for  the  PMMA  phantom  (compared  to  a  water 
phantom).  IMRT  plans  have  been  computed  using  a  water  phantom  and  the  PMMA  phantom  (all  have 
the  same  dimensions)  by  CORVUS  and  Monte  Carlo  and  the  results  were  compared  to  measurements. 
For  example,  for  a  9  field  IMRT  plan  with  15  MV  beams,  the  difference  between  the  calculated  and 
measured  dose  was  1.4%  for  the  water  phantom  and  1.6%  for  the  PMMA  phantom.  In  general,  the 
difference  between  the  measured  and  calculated  values  was  within  4%  with  an  estimated  measurement 
uncertainty  of  about  3%  and  Monte  Carlo  calculation  uncertainty  of  2%. 
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We  have  also  computed  dose  distributions  in  heterogeneous  PMMA  phantoms  with  either  lung  or  bone 
inserts.  Figure  5  shows  the  dose  distributions  calculated  by  Monte  Carlo  for  a  15  MV  beam  plan  and 
Table  I  shows  the  measured  and  calculated  dose  values  in  the  heterogeneous  phantom  for  both  4  and 
15  MV  beam  plans.  The  two  plans  were  generated  using  our  Monte  Carlo  based  inverse  planning 
system  using  the  same  parameters  except  for  the  beam  energy.  The  Monte  Carlo  calculated  values 
differed  from  the  measured  values  by  less  than  1%  for  both  4  and  15  MV  beam  plans  while  the 
CORVUS  values  differed  from  the  measured  values  by  1%  for  the  4  MV  plan  and  5%  for  the  15  MV 
plan.  Similar  agreement  was  obtained  for  other  phantom  configurations  and  beam  energies.  It  is 
concluded  that  the  Monte  Carlo  method  is  more  accurate  in  predicting  the  dose  in  heterogeneous 
phantoms  for  IMRT  dose  verification  and  can  be  used  to  validate  the  dose  distributions  for  IMRT  as 
part  of  the  QA  procedure^’. 


TABLE  I 


Dose  at  measurement  point  o  as  shown  in  Fig.  5. 


Beam  energy 

Monte  Carlo 

CORVUS 

Measurement 

4  MV 

2.177  Gy 

2.201  Gy 

2.177  Gy 

15  MV 

2.146  Gy 

2.276  Gy 

2.161  Gy 

5.  Investigate  beam  characteristics  of  electron  beams  modulated  by  magnetic  fields 

From  our  previous  studies,  it  seems  clear  that  an  electron  MLC  will  be  able  to  deliver  modulated 
electron  fields  accurately  and  therefore  He-filled  accelerator  design  may  not  be  necessary  (but  will 
further  improve  the  characteristics  if  implemented).  We  also  set  out  to  study  other  means  to  improve 
the  electron  dose  distributions  for  modulated  electron  beam  therapy  for  breast  cancer.  We  have 
achieved  some  preliminary  results  and  based  on  these  results  we  have  put  out  another  DOD  grant 
application  to  investigate  the  feasibility  of  using  magnetie  fields  to  modulate  the  beam  at  depth  to 
reduce  the  dose  to  the  distal  organs  such  as  the  lung  (see  below).  We  have  published  a  paper  on  this 
subject^*  and  we  briefly  describe  our  results  in  the  following  paragraphs. 

In  this  study,  Monte  Carlo  simulations  were  employed  to  study  the  characteristics  of  the  electron 
beams  of  a  clinical  linear  accelerator  in  the  presence  of  1.5  and  3.0  T  transverse  magnetic  fields  and  to 
assess  the  possibility  of  using  magnetic  fields  in  conjunction  with  modulated  eleetron  radiation  therapy 
(MERT).  The  starting  depth  of  the  magnetic  field  was  varied  over  several  centimeters.  It  was  found 
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that  peak  doses  of  as  much  as  2.7  times  the  surface  dose  could  be  achieved  with  a  1.5  T  magnetic 
field.  The  magnetic  field  was  shown  to  reduce  the  80%  and  20%  dose  drop-off  distance  by  50%  to 
80%.  The  distance  between  the  80%  dose  levels  of  the  pseudo-Bragg  peak  induced  by  the  magnetic 
field  was  found  to  be  extremely  narrow,  generally  less  than  1  cm.  However,  by  modulating  the  energy 
and  intensity  of  the  electron  fields  while  simultaneously  moving  the  magnetic  field,  a  homogeneous 
dose  distribution  with  low  surface  dose  and  a  sharp  dose  fall-off  was  generated.  Heterogeneities  are 
shown  to  change  the  effective  range  of  the  electron  beams,  but  not  eliminate  the  advantages  of  a  sharp 
depth-dose  drop-off  or  high  peak-to-surface  dose  ratio.  This  suggests  the  applicability  of  MERT  with 
magnetic  fields  in  heterogeneous  media.  The  results  of  this  study  demonstrate  the  ability  to  use 
magnetic  fields  in  MERT  to  produce  highly  desirable  dose  distributions. 

Monte  Carlo  calculations  were  performed  in  which  electrons  were  constrained  to  travel  down  the  same 
initial  axis,  i.e.  a  pencil  beam.  These  results  clearly  show  the  curvature  of  the  electron  beam  dose 
distribution  generally  along  the  expected  track.  By  superposition  of  many  pencil  beams  along  the 
lateral  axis,  it  is  clear  that  Bragg-peak-like  dose  peaks  can  be  constructed.  Narrow  fields,  defined  by  3 
X  3  cm  2  Cerrobend  cutouts  on  6  x  6  cm  2  applicators,  were  also  simulated.  What  is  apparent  is  that 
there  is  a  significant  shift  of  the  dose  peak  off  the  central  axis.  These  data  are  presented  primarily  as  a 
demonstration  of  the  effects  of  the  magnetic  field,  i.e.  the  curvature  of  the  beam  and  the  formation  of  a 
dose  peak.  However,  it  is  difficult  to  make  any  definitive  statements  regarding  depth-dose  distributions 
because  of  the  lack  of  equilibrium,  as  discussed  below. 

Previous  discussions  of  MERT  have  considered  the  use  of  fields  with  non-uniform  energies  and 
intensity  distributions.  Here,  the  technique  is  extended  to  the  use  of  multiple  magnetic  field  positions, 
i.e.  a  single  port  may  include  a  20  MeV  field  with  a  magnetic  field  beginning  at  3.0  cm  depth  for  a 
certain  number  of  monitor  units,  and  then  another  20  MeV  field  with  the  magnetic  field  located  at  a 
depth  of  5.0  cm  delivered  down  the  same  axis,  etc.  It  is  instructive  to  compare  this  methodology  with 
the  more  familiar  method  used  in  proton  beam  therapy.  In  that  system,  physical  blocks  of  different 
thicknesses  are  used  to  shift  the  Bragg  peak  to  different  depths,  and  intensity  modulation  is  provided 
by  allowing  the  beam  to  dwell  on  a  given  modulator  for  a  variable  time  period.  In  this  case,  the  only 
difference  is  the  means  of  moving  the  pseudo-Bragg  peak,  that  is,  by  moving  the  magnetic  field.  Beam 
intensity  is  modulated  simply  by  changing  the  number  of  monitor  units  delivered  with  the  magnetic 
field  at  a  given  position.  It  was  believed  that,  as  with  the  proton  beam,  useful  fields  can  be  constructed 
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while  maintaining  the  desirable  depth-dose  drop-off  and  low  skin  dose.  We  have  used  magnetic  fields 
modulated  electron  beamles  for  optimization  and  generate  treatment  plans  for  a  realistic  breast 
treatment  plan.  Figure  6  show  the  isodose  distributions  and  the  dose  volume  histogram  data  for  the  two 
plans,  one  with  magnetic  fields  and  the  other  without.  It  is  very  clear  that  magnetic  field  modulated 
electron  beams  show  much  better  target  conformity  and  normal  tissue  sparing.  Based  on  these  results, 
we  have  submitted  an  “IDEA”  proposal  to  further  investigate  this  novel  method,  and  if  funded  we  hope 
to  improvement  breast  MERT  significantly. 

We  can  demonstrate  how  this  works  as  follows.  A  simple  one-dimensional  optimization  routine  was 
utilized  to  generate  two  fields.  In  one  case,  a  target  profile  was  generated  where  the  dose  would  rise 
linearly  from  50%  at  the  surface  (relative  to  the  maximum  along  the  profiled  axis)  to  100%  at  1.0  cm, 
then  drop  to  zero  at  6.2  cm.  The  optimized  solution  was  solved  numerically,  and  a  linear  combination 
of  the  1.5  T  fields  was  generated  to  match  the  target  values  in  a  least-squares  sense.  The  results  are 
shown  in  figure  7.  Using  the  terminology  from  section  3.3,  the  width  of  the  treatment  region,  80/80, 
was  found  to  be  5.0  cm,  followed  by  a  80/20  drop-off  of  0.77  cm.  Between  the  depths  of  1.6  cm  and 
5.6  cm,  the  dose  was  homogeneous  to  within  5  4%  of  the  average  dose.  For  more  detailed  descriptions 
of  the  magnetic  field  modulation  please  see  reference  [28]. 

Key  Research  Accomplishments 

We  have  accomplished  the  following  tasks: 

•  Develop  fast  iterative  optimization  algorithms:  A  fast  iterative  optimization  algorithm  has  been 
developed  for  electron  beam  optimization.  The  algorithm  has  been  tested  for  both  electron 
treatment  planning  optimization  and  photon  beam  optimization. 

•  Install  and  test  different  objective  function  for  electron  beam  optimization:  We  have  installed 
several  objective  functions  for  the  optimizer  used  for  electron  beam  EEMRT.  We  have  selected  a 
dose  based  objective  function  with  soft  constraints  using  dose  volume  histogram  information. 

•  Study  the  characteristics  of  electron  beam  modulation  using  magnetic  fields:  We  have  studied  the 
characteristics  of  small  field  electron  beams  modulated  by  an  electron  MLC  and  0.5-3.0  T 
magnetic  fields.  We  have  studied  beam  optimization  using  these  modulated  electron  beamlets. 
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•  Verify  electron  beam  ElMRT  using  EGS4/MCD0SE:  We  have  developed  an  EGS4  user  code  called 
MCDOSE  to  calculate  dose  distribution  for  electron  treatment  planning.  We  have  compared  the 
dose  distributions  calculated  using  MCDOSE  and  those  by  DOSXYZ  and  achieved  good 
agreement  (within  1%  of  maximum  dose). 

•  Verify  photon  IMRT  using  EGS4/MCDOSE:  We  have  used  MCDOSE  to  calculate  dose  distribution 
for  photon  IMRT  treatment  planning.  We  have  compared  the  dose  distributions  calculated  using 
MCDOSE  and  measurements  and  achieved  good  agreement  (within  3%  of  maximum  dose). 

Reportable  Outcomes 

Peer-reviewed  papers  resulting  from  or  supported  in  part  by  this  grant: 

•  C.M.  Ma,  T.  Pawlicki,  M.C.  Lee,  S.B.  Jiang,  J.S.  Li,  J.  Deng,  E.  Mok,  B.  Yi,  G.  Luxton  and  A.L. 
Boyer,  Energy-  and  intensity-modulated  electron  beams  for  radiotherapy,  Phys.  Med.  Biol.  (2000) 
45:2293-2311 

•  J.S.  Li,  T.  Pawlicki,  J.  Deng,  S.B.  Jiang,  E.  Mok  and  C.M.  Ma,  Validation  of  a  Monte  Carlo  dose 
calculation  tool  for  radiotherapy  treatment  planning,  Phys.  Med.  Biol.  (2000)  45:  2969-2985 

•  S.B.Jiang,  A.  Kapur  and  C.M.  Ma,  Electron  beam  modelling  and  commissioning  for  Monte  Carlo 
treatment  planning,  Med.  Phys.  (2000)  27:180-191 

•  C.M.  Ma,  T.  Pawlicki,  S.B.  Jiang,  E.Mok,  A.  Kapur,  L.  Xing,  L.  Ma  and  A.L.  Boyer,  Monte  Carlo 
verification  of  IMRT  dose  distributions  from  a  commercial  treatment  planning  optimization 
system,  Phys.  Med.  Biol.,  (2000)  45:2483-2495 

•  M.  C.  Lee  and  C.M.  Ma,  Monte  Carlo  charactarization  of  clinical  electron  beams  in  transverse 
magnetic  fields,  Phys.  Med.  Biol.  (2000)  45:2947-2967 

Meeting  abstracts  resulting  from  or  supported  in  part  by  this  grant: 

•  C.M.  Ma,  JS  Li,  T  Pawlicki,  SB  Jiang  and  J  Deng,  MCDOSE  -  a  Monte  Carlo  dose  calculation  tool 
for  radiation  therapy  treatment  planning,  Proc.  of  the  XIII  International  Conference  on  the  Use  of 
Computer  in  Radiation  Therapy  (ICCR),  Eds:  W  Schlegel  and  T  Bortfeld  (Springer- Verlag, 
Heidelberg)  2000.  123-125 

•  C.M.  Ma,  T  Pawlicki,  MC  Lee,  SB  Jiang,  JS  Li,  J  Deng,  E  Mok  and  AL  Boyer,  Modulated  electron 
beams  for  treatment  of  breast  cancer.  Oral  presentation  at  the  XIII  International  Conference  on  the 
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Use  of  Computer  in  Radiation  Therapy  (ICCR),  Heidelberg,  Germany,  May  22  -  25,  2000. 
Proceedings  to  be  published  at  the  meeting. 

•  C.M.  Ma,  Monte  Carlo  methods  in  electron  beams  treatment  planning,  invited  talk  at  the  41st 
Annual  Meeting  of  the  American  Association  of  Physicists  in  Medicine,  Nashville,  TN,  July  24-29, 
1999.  Med.  Phys.,  26(1999  AAPM  Annual  Meeting  Program)  1999. 

•  J.  Deng,  S.B.  Jiang,  T.  Pawlicki,  J.  Li  and  C.M.  Ma,  Electron  beam  commissioning  for  Monte 
Carlo  dose  calculation,  ICCR  2000. 

•  C.M.  Ma,  T  Pawlicki,  SB  Jiang,  JS  Li,  Deng,  D  Findley,  E  Mok  and  AL  Boyer,  Implementation  of 
a  Monte  Carlo  dose  calculation  module  in  the  FOCUS  treatment  planning  system,  CMS’ FOCUS 
2000  User’s  Symposium,  St.  Louis,  MO,  April  10-11, 2000. 

•  J.  S.  Li,  T.  Pawlicki,  J.  Deng,  S.B.  Jiang,  A.  Kapur,  E.  Mok  and  C.M.  Ma,  Clinical 
validation  of  a  Monte  Carlo  dose  calculation  code  for  radiotherapy  treatment  planning 
(abstract),  Med.  Phys.  26,  1083 

•  J.  S.  Li,  T.  Pawlicki,  J.  Deng,  S.B.Jiang  and  C.M.  Ma,  Simulation  of  Beam  Modifiers  for 
Monte  Carlo  Treatment  Planing,  submitted  to  ICCR  2000 

•  S.  B.  Jiang,  A.  Kapur,  and  C.M.  Ma,  Electron  beam  modeling  and  eommissioning  for 
Monte  Carlo  treatment  planning,  Med.  Phys.  26,  1084  (1999)  (abstract),  AAPM  Annual 
Meeting  (Nashville). 

Funding  applied  for  based  on  work  resulting  from  or  supported  in  part  by  this  grant: 

1.  US  Army  Breast  Cancer  Research  Program:  Postdoctoral  Training  Grant  (PI:  S.B.  Jiang;  Mentor: 
C.-M.  Ma):  Investigation  of  an  electron  MLC  for  energy  and  intensity  modulated  electron  beam 
radiotherapy  of  breast  cancer  (awarded  in  2000) 


2.  US  Army  Breast  Cancer  Research  Program:  Postdoctoral  Training  Grant  (PI:  T.  Pawlicky;  Mentor: 
C.-M.  Ma):  Effect  of  organ  motion  and  patient  setup  uncertainty  on  the  treatment  of  breast  cancer 
with  energy  and  intensity  modulated  electron  beam  radiotherapy  (awarded  in  2000) 


3.  US  Army  Breast  Cancer  Research  Program:  Predoctoral  Training  Grant  (PI:  Michael  C.  Lee; 

Mentor;  C.-M.  Ma):  Beam  verification  for  modulated  electron  radiaiton  therapy  treatment  of  breast 
cancer  (submitted  in  2(X)0) 


4.  US  Army  Breast  Cancer  Research  Program:  Idea  Grant  (PI:  C.-M.  Ma):  Electron  beam  modulation 
using  magnetic  fields  for  breast  cancer  treatment  (submitted  in  2000) 
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Conclusions 


We  have  made  significant  progress  during  our  second-year  investigation.  We  have  successfully 
performed  the  tasks  scheduled  in  the  “Statement  of  Work”.  We  have  developed  optimization 
algorithms  for  comparisons  of  the  optimized  treatment  plans  with  photon  IMRT  or  electron  EIMRT. 
We  have  installed  and  investigated  suitable  objective  functions  for  electron  beam  optimization.  We 
have  verified  electron  EIMRT  and  photon  IMRT  using  the  EGS4/MCDOSE  system  and  compared 
dose  distributions  for  conventional  treatment  modalities.  We  also  studied  the  effect  of  magnetic  fields 
on  modulation  of  electron  beam  profile  and  depth  dose  curves  for  EIMRT  of  breast  cancer. 

“So  what?” 

Our  second  year  results  have  provided  evidence  to  support  the  hypothesis  of  this  proposal  that  by  using 
computer-controlled,  specially  designed  multi-leaf  collimators  for  electron  beams  and  the  Monte  Carlo 
treatment  planning  technique,  energy-  and  intensity-modulated  electron  beams  may  be  optimized  to 
significantly  improve  the  dose  uniformity  in  the  target  volume,  to  reduce  the  dose  to  the  critical 
structures  nearby  and  therefore  reduce  the  late  effects  associated  with  breast  cancer  radiotherapy. 
Further  studies  as  scheduled  in  the  “Statement  of  Work”  are  needed  to  investigate  treatment  plans  for 
realistic  breast  cases  generated  using  different  beam  modalities.  We  shall  be  able  to  finally  determine 
the  degree  of  feasibility  of  using  electron  EIMRT  to  improve  breast  cancer  treatment. 
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Appendices 

List  of  Figures  quoted  in  the  body  of  text: 


Figure  1  Comparisons  of  the  lateral  dose  profiles  at  different  depths  in  a  water  phantom  between 
MCDOSE  and  EGS4/BEAM  plus  DOSXYZ.  A  45°  wedge  were  simulated  for  a  15MV  photon  beam, 
the  field  size  was  10cm  x  10cm  defined  at  lOOSSD. 
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Figure  2  Depth  dose  curves  (a)  in  a  Tissue-Lung-Tissue  phantom  calculated  by  DOSXYZ  and 
MCDOSE  for  6MeV,  15MV  and  2MeV  photon  beams  with  field  size  6  cm  x  6  cm  defined  at  lOOSSD 
and  lateral  dose  profiles  at  depths  of  3.3cm,  9.0cm  and  15.0cm  for  15MV  photon  beam  (b).  From  depth 
of  9.0cm  to  12.5cm,  it  was  material  of  lung  from  -5.0cm  to  5.0cm  in  X  dimension  and  from  -2.0cm  to 
2.0cm  in  Y  dimension. 


(a) 


X  Axis  (cm  ) 


(b) 
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Fig.  3  Depth  dose  curves  (a)  in  a  Tissue-Bone-Tissue  phantom  calculated  by  DOSXYZ  and  MCDOSE 
for  6MeV,  15MV  and  2MeV  photon  beams  with  field  size  6  cm  □  6  cm  defined  at  lOOSSD  and  lateral 
dose  profiles  at  depths  of  3.3cm,  8.0cm  and  11.0cm  for  15MV  photon  beam  (b).  From  depth  of  7.0cm 
to  10.0cm,  it  was  material  of  bone  from  -2.0cm  to  2.0cm  in  Y  dimension  and  from  -5.0cm  to  5.0cm  in 
X  dimension. 


(a) 


Y  Axis  (cm  ) 


(b) 
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Figure  4  The  beam  eye  views  of  the  dose  distributions  at  DMAX  in  a  water  phantom  for  a  10cm  x 
10cm  12MeV  electron  beam.  The  beam  was  modified  using  an  electron  cutout  with  a  butterfly-shaped 
opening.  The  simulation  results  of  BEAM  (thin  lines)  and  MCDOSE  (thick  lines)  are  shown. 


Figure  5  Dose  distribution  for  a  hypothetical  target  in  the  PMMA  phantom  with  bone  insert  calculated 
by  Monte  Carlo.  Nine  15  MV  photon  beams  were  used.  The  dose  at  the  center  of  the  target  was 
measured  using  an  ion  chamber. 


Figure  6  (a)  Isodose  curves  for  a  breast  treatment  plan.  Thick  lines  are  optimized  dose  distributions  for 
MERT  with  1.5  T  magnetic  fields  and  thin  lines  are  MERT  without  magnetic  fields.  The  target 
coverage  is  similar  but  the  dose  to  the  lung  behind  the  target  is  almost  completely  removed  with  the 
“soft  collimation”  using  a  magnetic  field,  (b)  dose  volume  histograms  for  the  two  plans  shown  in  (a). 


(a) 

DVH  Data 


(b) 
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Figure  7  Central-axis  depth-dose  curves  demonstrating  how  linear  combinations  of  10  x  10  cm  2 
electron  beams  in  magnetic  fields,  with  different  intensities  and  magnetic  field  placements,  can  be  used 
to  create  homogeneous  dose  distributions.  The  dose  curves  were  for  the  1.5  T  optimization. 
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Abstract.  This  work  investigates  the  feasibility  of  optimizing  energy-  and  intensity-modulated 
electron  beams  for  radiation  therapy.  A  multileaf  collimator  (MLC)  specially  designed  for 
modulated  electron  radiotherapy  (MERT)  was  investigated  both  experimentally  and  by  Monte  Carlo 
simulations.  An  inverse-planning  system  based  on  Monte  Carlo  dose  calculations  was  developed  to 
optimize  electron  beam  energy  and  intensity  to  achieve  dose  conformity  for  target  volumes  near  the 
surface.  The  results  showed  that  an  MLC  with  5  mm  leaf  widths  could  produce  complex  field  shapes 
for  MERT.  Electron  intra-  and  inter-leaf  leakage  had  negligible  effects  on  the  dose  distributions 
delivered  with  the  MLC,  even  at  shallow  depths.  Focused  leaf  ends  reduced  the  electron  scattering 
contributions  to  the  dose  compared  with  straight  leaf  ends.  As  anticipated,  moving  the  MLC 
position  toward  the  patient  surface  reduced  the  penumbra  significantly.  There  were  significant 
differences  in  the  beamlet  distributions  calculated  by  an  analytic  3-D  pencil  beam  algorithm  and 
the  Monte  Carlo  method.  The  Monte  Carlo  calculated  beamlet  distributions  were  essential  to  the 
accuracy  of  the  MERT  dose  distribution  in  cases  involving  large  air  gaps,  oblique  incidence  and 
heterogeneous  treatment  targets  (at  the  tissue-bone  and  bone-lung  interfaces).  To  demonstrate  the 
potential  of  MERT  for  target  dose  coverage  and  normal  tissue  sparing  for  treatment  of  superficial 
targets,  treatment  plans  for  a  hypothetical  treatment  were  compared  using  photon  beams  and  MERT. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version;  see  www.iop.org) 


1.  Introduction 

Photon  beams  have  been  an  effective  modality  for  breast  cancer  treatment  in  radiation  therapy. 
Although  such  conventional  treatment  with  tangential  photon  fields  has  been  successful,  the 
following  two  problems  (or  potential  areas  of  improvement)  remain: 

(a)  The  inclusion  of  the  lung  and  other  normal  tissues,  and  sometimes  of  a  small  volume  of 
the  heart,  in  the  high-dose  volume  due  to  tumour  location,  patient  size  or  in  the  case  of 
chest-wall  treatments. 

(b)  High  exit  or  scatter  dose  to  the  normal  structures  such  as  the  lung,  the  heart  and  the 
contralateral  breast. 

Advances  in  the  state  of  the  art  of  computer-controlled  medical  linear  accelerators  have 
recently  become  available  that,  along  with  newly  developed  treatment  planning  techniques, 
may  provide  significant  improvements  in  the  delivery  and  control  of  external  beam  radiation 
through  beam-intensity  modulation  (Boesecke  et  al  1988,  Brahme  1988,  Convery  and 
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Rosenbloom  1992,  Leibel  et  al  1992,  Webb  1992,  1997,  LoSasso  et  al  1993,  Powlis  et  al 

1993,  Chui  et  al  1994,  Mageras  et  al  1994,  Brewster  et  al  1995,  Fraass  et  al  1995,  Kutcher 
et  al  1995,  Mackie  et  al  1995,  McShan  et  al  1995,  Ling  et  al  1996,  Boyer  et  al  1997).  It  is 
expected  that  using  photon  IMRT,  the  problem  (a)  above  may  be  significantly  improved  but  (b) 
may  become  more  serious  as  treatment  time  increases  with  the  number  of  fields/segments  used 
(increased  leakage  or  scattering  dose).  Using  the  modulated  electron  radiotherapy  (MERT) 
technique  (Lief  etal  1996,  Hyodynmaa  et  al  1996,  Zackrisson  and  Karlsson  1996,  Asell  etal 
1997,  Ebort  and  Hoban  1997,  Karlsson  etal  1998,  1999),  on  the  other  hand,  problem  (a)  may 
also  be  significantly  improved  and  problem  (b)  may  almost  be  eliminated  due  to  the  nature  of 
the  electron  beams. 

In  the  optimization  process  of  MERT,  dose  conformity  along  the  beam  direction  can  be 
achieved  by  modulating  the  electron  incident  energy,  making  use  of  the  sharp  dose  fall-off 
feature.  A  drawback  is  its  large  penumbra  at  large  depths.  Traditionally,  electron  beams  are 
shaped  using  a  cutout  (or  blocks)  and  beam  penetration/intensity  may  be  modified  using  a  bolus. 
However,  it  is  time-consuming  to  make  such  beam  modifiers  and  the  treatment  time  would  be 
significantly  increased  if  such  beam  modifiers  are  used  for  MERT.  Efforts  have  been  made  to 
use  computer-controlled  MLC  for  electron  beam  modulation.  The  recent  results  by  Karlsson 
et  al  (1999)  showed  that  by  replacing  the  air  in  the  treatment  head  with  a  low-cost,  custom- 
made  helium  balloon,  the  beam  penumbral  width  (20/80)  was  reduced  from  18  to  11  mm  at 
80  cm  SSD.  The  beam  characteristics  are  affected  by  the  position  of  the  MLC.  However,  by 
replacing  the  air  between  the  MLC  and  the  patient  with  a  helium  balloon,  the  beam  penumbra 
become  almost  the  same  as  that  achieved  by  electron  beam-shaping  with  an  electron  applicator 
that  extends  to  the  patient  skin  surface  (Karlsson  et  al  1999).  This  means  that  many  of  the 
techniques  so  far  developed  with  computer-controlled  MLC  and  our  experience  with  MLC 
photon  beam  modulation  can  be  adopted  for  use  with  MERT. 

The  calculation  of  dose  distributions  for  electron  beam  radiotherapy  planning  is 
challenging  because  electron  scattering  is  strongly  affected  by  changes  in  density  and 
composition  in  the  patients.  The  3D  pencil  beam  algorithm  (Hogstrom  et  al  1981)  is  a  fast 
analytical  algorithm  which  has  been  adopted  by  most  treatment  planning  systems.  However, 
it  has  limitations  with  small  irregular  electron  fields  and  in  the  presence  of  inhomogeneities 
(Cygler  et  al  1987,  Bielajew  et  al  1987,  Mah  et  al  1989,  Mackie  et  al  1994,  Ma  et  al  1999). 
The  Monte  Carlo  simulation  has  been  demonstrated  to  be  a  viable  option  for  such  complex 
situations,  and  also  the  only  way  to  take  into  account  back-scattering  from  denser  materials  in 
a  patient  (e.g.  bone  or  metal  inserts)  (Shortt  etal  1986,  Cygler  etal  1987,  Mackie  etal  1994, 
Kawrokaw  etal  1996,  Mohan  1997,  Kapur  1999,  Ma  et  al  1999).  The  EGS4/BEAM  system 
was  developed  for  the  simulation  of  radiotherapy  beams  from  various  radiotherapy  treatment 
units,  such  medical  accelerators  (Rogers  et  al  1995).  Excellent  agreement  (1-3%)  has  been 
achieved  between  the  Monte  Carlo  dose  distributions  calculated  using  the  simulated  particle 
phase-space  data  and  measurements  (Rogers  et  al  1995,  Kapur  et  al  1998,  Zhang  et  al  1999, 
Ma  et  al  1999).  We  have  installed  a  Monte  Carlo  patient  dose  calculation  tool  on  a  clinical 
treatment  planning  system  (Ma  etal  1999)  and  used  Monte  Carlo  for  treatment  planning  and 
dose  delivery  validation.  This  has  reduced  the  uncertainty  of  the  accelerator  output  for  small 
irregular  field  electron  beams  from  up  to  10%  to  about  3%  (Ma  et  al  1997,  Kapur  et  al  1998). 

Conformal  radiotherapy  was  initially  used  to  limit  the  normal  tissue  dose  by  conforming 
the  treatment  field  to  the  beam’s-eye-view  projection  of  the  target  volume  (Takahashi  1965). 
For  photon  beams,  the  MLC  was  used  to  collimate  the  fields  and  later  to  modulate  the  beam 
intensity  in  the  field  (Boesecke  et  al  1988,  Brahme  1988,  Convery  and  Rosenbloom  1992, 
Leibel  et  al  1992,  Webb  1992,  1997,  LoSasso  et  al  1993,  Powlis  et  al  1993,  Chui  et  al 
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Figure  1.  A  prototype  of  an  electron  MLC  mounted  on  the  bottom  scraper  of  a  25  cm  x  25  cm 
applicator  on  a  Varian  Clinac  2100C  accelerator.  The  MLC  has  30  pairs  of  steel  leaves  and  the  leaf 
positions  are  fixed  by  the  tightening  screws. 


McShm  etal  1995,  Ling  et  al  1996,  Boyer  et  al  1997).  There  have  been  a  few  studies 
on  electron  conformal  therapy.  Tailoring  dose  distributions  using  electron  beams  requires 
substantial  beam  manipulation,  due  to  their  scattering  characteristics.  Such  manipulation  is 
already  possible  with  radiation  sources  such  as  microtrons  where  preferential  energy  selection 
and  magnetically  scanned  pencil  beams  are  possible  (Lief  et  al  1996).  Both  intensity-  and 
energy/intensity-modulated  electron  beams  have  been  investigated  to  conform  the  dose  to  the 
target  near  the  surface  (Lief  etal  1996,  Hyodynmaa  a/ 1996,  Zackrisson  and  Karlsson  1996, 
Ebert  and  Hoban  1997,  Karlsson  etal  1998).  More  recent  work  has  studied  the  combination  of 
photon  IMRT  and  MERT  for  targets  at  greater  depth  (Karlsson  etal  1999).  Using  the  helium- 
balloon  technique  together  with  a  computer-controlled  MLC,  it  may  be  possible  to  deliver  a 
set  of  intensity-modulated  beams  with  different  energies  and  incident  angles. 

In  this  work,  we  investigate  the  feasibility  of  optimizing  energy-  and  intensity-modulated 
electron  beams  for  radiotherapy  treatment.  We  report  here  our  Monte  Carlo  studies  of 
a  multileaf  collimator  specially  designed  for  MERT  and  some  preliminary  experimental 
results.  We  also  report  on  the  dose  calculation  algorithms  and  their  effects  on  treatment  plan 
optimization  for  MERT.  We  will  discuss  the  differences  in  the  beam  characteristics  between 
a  photon  MLC  and  an  electron  MLC.  We  will  compare  the  dose  distributions  between  a 
conventional  tangential  photon  treatment  plan  and  a  MERT  treatment  plan  for  a  hypothetical 
breast  treatment  to  demonstrate  the  potential  of  MERT  for  target  dose  coverage  and  normal 
tissue  sparing. 

2.  Materials  and  method 

2.7.  The  prototype  electron  MLC 

We  have  developed  a  prototype  manual-driven  electron  MLC  for  the  beam  delivery  for  MERT. 
As  shown  in  figure  1,  the  electron  MLC  consists  of  30  steel  leaf  pairs,  which  were  made  from 
the  off-the-shelf  steel  bars  for  convenience  and  cost-effectiveness.  Each  leaf  is  0.476  cm  wide, 
20  cm  long  and  2.54  cm  thick  with  straight  edges  and  ends.  The  leaves  were  mounted  on  a 
steel  frame,  which  can  be  attached  to  the  bottom  scraper  of  a  25  cm  x  25  cm  electron  applicator 
on  a  Varian  Clinac  2100C  (Varian  Oncology  Systems,  Palo  Alto,  CA).  The  leaves  can  slide  in 
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the  steel  frame  and  the  leaf  positions  can  be  easily  set  using  a  pre-cut  cardboard  for  a  beam 
segment.  The  field  shape  is  maintained  by  tightening  the  screws  from  the  side.  The  largest 
radiation  field  available  using  the  electron  MLC  was  15.7  cm  x  15.7  cm  projected  at  100  cm 
source-surface  distance  (SSD). 

Because  of  the  existing  electronic  device  for  inserting  the  electron  cutout,  the  leaves  could 
not  be  placed  at  the  last  scraper  level  without  modifying  the  existing  applicator  geometry. 
Instead,  the  steel  frame  was  inserted  using  the  electron  cutout  mount  and  the  leaves  were 
placed  immediately  above  the  last  scraper.  This  resulted  in  a  slightly  greater  air  gap  (7  cm) 
between  the  bottom  of  the  electron  MLC  leaves  and  the  phantom  surface  (assuming  a  97  cm 
SSD)  compared  with  that  of  an  electron  cutout  (5  cm  for  a  100  cm  SSD).  The  projected  leaf 
width  for  a  5  cm  air  gap  is  0.5  cm,  while  for  a  7  cm  air  gap  the  projected  leaf  width  is  0.51  cm 
(e.g.  for  the  current  configuration  at  97  cm  SSD).  Further  modifications  are  needed  to  the 
electron  applicator  geometry  in  order  to  lower  the  electron  MLC  leaves.  The  ideal  location 
for  the  MLC  leaves  is  the  last  scraper  since  electron  cutouts  will  no  longer  be  needed  if  an 
electron  MLC  is  in  place. 

Film  measurement  was  performed  to  study  the  characteristics  of  the  electron  beams 
collimated  by  the  electron  MLC.  The  film  was  calibrated  following  the  AAPM  TG-25 
recommendations  (AAPM  1991)  and  the  exposures  were  taken  by  placing  film  at  different 
depths  in  a  solid  water  phantom.  The  film  was  scanned  using  a  film  scanner  which  has  a 
spatial  resolution  of  about  0.15  cm. 

2.2.  The  Monte  Carlo  beam  simulation 

We  have  used  the  EGS4  (Nelson  et  al  1985)  user  code  BEAM  for  the  accelerator  head 
simulation.  Detailed  descriptions  of  the  software  can  be  found  in  Rogers  et  al  (1995).  A 
detailed  description  of  the  clinical  implementation  of  the  Monte  Carlo  method  at  the  Stanford 
Medical  Center  was  given  in  a  previous  publication  (Ma  et  al  1999). 

For  this  work,  we  have  used  the  previously  simulated  Monte  Carlo  beam  data  for  6,  12 
and  20  MeV  electron  beams  from  a  Varian  Clinac  2100C  linear  accelerator  and  for  6  MV 
photon  beams  from  a  2300CD  accelerator  (Varian  Oncology  Systems,  Palo  Alto,  CA).  The 
dimensions  and  materials  for  the  accelerator  components  were  incorporated  according  to  the 
manufacturer’s  specifications.  Electron  beams  emerging  from  the  vacuum  exit  window  were 
assumed  to  be  monoenergetic  and  monodirecti onal  with  a  beam  radius  of  0. 1-0.2  cm  (Kapur 
et  al  1998).  The  energy  cutoffs  for  electron  transport  in  the  accelerator  simulation  (ECUT 
and  AE)  were  700  keV  (kinetic  +  rest  mass)  and  for  photon  transport  (PCUT  and  AP)  10  keV. 
The  electron  transport  step  length  was  confined  such  that  the  maximum  fractional  energy  loss 
per  electron  step  is  4%  (i.e.  ESTEPE  =  0.04).  The  ICRU  recommended  compositions  and 
stopping  power  values  were  used  for  the  materials  in  the  accelerator  simulations  (ICRU  1984). 
The  phase-space  data  were  scored  at  a  plane  either  immediately  above  the  photon  MLC  or 
above  the  lowest  scraper.  The  number  of  particles  was  about  2-30  million  in  an  electron  beam 
file  and  about  50  million  in  a  photon  file. 

Field  shaping  by  the  photon  MLC  or  electron  MLC  was  further  simulated  using  the  BEAM 
component  module  MLC.  MLC  could  simulate  either  straight  or  ‘double  focused’  leaf  edges 
and  ends.  In  this  work,  we  have  simulated  electron  beams  collimated  by  a  photon  MLC  with 
both  straight  and  double  focused  MLC  leaf  shapes.  The  leaves  were  7.5  cm  thick  and  made  of 
tungsten.  The  leaf  center  was  49  cm  from  the  isocentre.  The  intervening  air  in  the  accelerator 
and  between  the  MLC  and  the  isocentre  was  in  some  cases  replaced  with  helium  to  investigate 
the  effect  of  electron  scattering  in  the  air.  In  the  simulations  of  the  electron  beams  collimated  by 
an  electron  MLC,  the  leaves  were  placed  on  the  bottom  scraper  of  a  25  cm  x  25  cm  applicator 
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with  a  7  cm  air  gap  between  the  bottom  surface  of  the  MLC  and  the  isocentre.  Tungsten  leaves 
of  1.5  cm  thickness  with  straight  edges  and  ends  were  used  in  all  the  simulations  and  the  phase 
space  data  were  used  in  the  subsequent  dose  calculations  except  for  the  leaf  leakage  study 
where  different  leaf  materials  and  thicknesses  were  investigated  for  the  electron  MLC. 


23.  The  Monte  Carlo  dose  calculation 

The  EGS4  user  code,  MCDOSE  (Ma  etal  1 999),  was  used  in  this  work  for  the  dose  calculations. 
MCDOSE  was  designed  for  dose  calculations  in  a  3D  rectilinear  voxel  geometry.  Voxel 
dimensions  were  completely  variable  in  all  three  directions.  Every  voxel  (volume  element) 
could  be  assigned  to  a  different  material.  The  cross-section  data  for  the  materials  used  were 
available  in  a  pre-processed  PEGS4  cross-section  data  file.  The  mass  density  of  the  material 
in  a  MCDOSE  calculation  was  varied  based  on  the  patient’s  CT  data  although  the  density 
effect  corrections  for  the  stopping  powers  of  the  material  remain  unchanged  (Ma  et  al  1999). 
The  voxel  dimensions  and  materials  were  defined  in  a  MCDOSE  input  file  together  with  the 
transport  parameters  such  as  the  energy  cutoffs  (ECUT  and  PCUT),  the  maximum  fractional 
energy  loss  per  electron  step  (ESTEPE),  and  the  parameters  required  by  PRESTA  (Bielajew 
and  Rogers  1987).  Several  variance  reduction  techniques  have  been  implemented  in  the 
MCDOSE  code  to  improve  the  calculation  efficiency.  These  include  photon  interaction  forcing, 
particle  splitting,  Russian  roulette,  electron  range  rejection  and  region  rejection,  particle  track 
displacement  and  rotation,  and  correlated  sampling.  Detailed  descriptions  of  these  techniques 
have  been  given  elsewhere  (Rogers  and  Bielajew  1990,  Ma  and  Nahum  1993,  Rogers  et  al 
1995,  Kawrakow  et  al  1996,  Keall  and  Hoban  1996,  Ma  et  al  1999). 

For  patient  dose  calculations,  the  simulation  phantom  was  built  from  the  patient’s  CT  data 
with  up  to  128  X  128  X  128  voxels  (uniform  in  any  dimension).  The  side  of  a  voxel  varied 
from  0.2  to  0.4  cm.  A  separate  program  was  developed  to  convert  the  patient’s  CT  data  from 
the  FOCUS  treatment  planning  system  (Computerized  Medical  Systems,  St  Louis,  MO)  to 
desired  dimensions,  material  types  and  densities.  The  organ  contours  were  also  obtained  for 
dose  calculation  and  analysis.  The  phase-space  data  obtained  from  a  BEAM  simulation  were 
used  as  a  source  input  with  variable  source  positions  and  beam  incident  angles.  To  simulate 
the  dose  distribution  of  a  finite  size  beamlet  used  by  the  inverse  planning  process,  particles 
were  transported  to  the  MLC  plane  and  only  those  within  the  beamlet  area  (=  1  cm  x  1  cm 
projected  at  100  cm  SSD)  were  allowed  to  go  through.  This  ignored  the  bremsstrahlung 
photon  leakage  and  electron  scattering  by  the  leaf  ends  in  the  optimization  process  (the  effect 
was  corrected  in  the  final  dose  calculation,  as  discussed  below).  After  optimization,  a  leaf 
sequence  was  generated  using  a  modified  ‘step  and  shoot’  algorithm  based  on  our  early  work 
(Ma  et  al  1998).  The  final  MERT  dose  distribution  was  computed  based  on  an  intensity  map 
(a  2D  distribution  of  particle  weighting  factors)  reconstructed  from  the  leaf  sequence.  The 
bremsstrahlung  leaf  leakage  effect  was  included  in  the  intensity  map  using  the  leaf  sequence 
and  pre-calculated  leaf  leakage  data  for  1.5  cm  thick  tungsten  leaves.  MCDOSE  produced 
data  files  that  contained  geometry  specifications  such  as  the  number  of  voxels  in  all  the  three 
directions  and  their  boundaries  as  well  as  the  dose  values  and  the  associated  (la)  statistical 
uncertainties  in  the  individual  voxels  and  organs  (structures).  The  EGS4  transport  parameters 
were  ECUT  =  AE  =  700  keV,  PCUT  =  AP  =  10  keV  and  ESTEPE  =  0.04.  The  number  of 
particle  histories  simulated  ranged  from  2  million  to  30  million  for  a  MERT  treatment.  The 
lor  statistical  uncertainty  in  the  dose  was  generally  2%  or  smaller  of  the  Dmux  value.  The  CPU 
time  required  for  a  MERT  simulation  was  about  1-3  h  on  a  Pentium  III  450  MHz  PC  with  the 
variance  reduction  option  switched  on. 
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Figure  2.  Beam  intensity  distributions  measured  by  film  on  the  surface  of  a  solid  water  phantom 
for  6,  12  and  20  MeV  electrons.  The  MLC  leaf  positions  for  the  electron  fields  are  also  shown 
(bottom  left). 


2.4.  The  optimization  process 

The  treatment  planning  optimization  system  used  in  this  work  is  a  home-developed  system 
based  on  the  work  by  Jiang  (1998).  First,  the  planner  inputs  the  patient  geometry  and  defines 
the  treatment  setup,  such  as  the  beam  energy,  number  and  orientations  of  beams,  etc.  The 
target  volume  and  the  critical  structures  are  defined  by  the  clinician.  A  reference  monitor  unit  is 
assigned  to  each  open  rectangular  beam  and  the  dose  deposition  coefficients,  which  are  defined 
as  the  dose  contribution  from  a  beamlet  to  a  point,  are  calculated  using  the  MCDOSE  code. 

Second,  using  the  calculated  dose  deposition  coefficients  as  input,  the  optimal  intensity 
profile  for  each  beam  is  achieved  using  a  gradient  method  to  minimize  the  objective  function. 
For  the  target  area,  a  quadratic  form  of  objective  function  is  specified.  In  addition,  two 
target  dose-uniformity  constraints  are  used  to  ensure  a  uniform  target  dose  distribution  and 
to  distinguish  the  clinical  importance  of  cold  and  hot  spots.  For  the  critical  structures, 
maximum-dose  constraint  and  several  levels  of  dose-volume  constraints  are  assigned  to  each 
structure.  For  each  objective  function  and  constraint,  an  importance  weight  relative  to  the 
target  objective  function  is  assigned.  All  the  constraints  are  mathematically  transformed  to 
the  penalty  functions  of  quadratic  forms.  The  augmented  objective  function,  which  should  be 
minimized,  is  a  combination  of  the  original  objective  functions  and  all  penalty  functions.  The 
results  of  the  optimization  process  are  the  intensity  profiles  for  the  individual  fields  (different 
incident  energies  and  gantry  angles).  The  same  optimizer  has  been  used  for  photon  beam 
optimization  with  the  Monte  Carlo  method  and  a  finite-size  pencil  beam  algorithm  (Jiang 
1998,  Jiang  etal  1999,  Pawlicki  etal  1999). 

3.  Resultsand  discussion 

3.1.  Characteristics  of  electron  beams  collimated  by  an  electron  MLC 

Figure  2  shows  the  electron  fields  collimated  by  the  prototype  electron  MLC  for  6,  12  and 
20  MeV  electron  beams  on  a  Varian  Clinac  2100C  machine.  For  convenience,  a  photo  showing 
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Figures.  Beam  profiles  measured  by  film  on  the  surface  of  a  solid  water  phantom  for  the  6  MeV 
and  20  MeV  electron  fields  shown  in  figure  2:  (a)  20  MeV  along  A-A;  (b)  20  MeV  along  B-B; 
(c)  6  MeV  along  A-A;  (d)  6  MeV  along  B-B. 


the  MLC  leaf  positions  for  the  field  shape  is  also  included  in  figure  2.  Figure  2  shows  the 
film  measurement  at  the  surface  of  a  solid  water  phantom  (97  cm  SSD)  for  6,  12  and  20  MeV 
electron  beams.  Figure  3  shows  the  measured  profiles  on  the  phantom  surface  along  A-A  and 
B-B  for  the  6  MeV  and  20  MeV  electron  fields  shown  in  figure  2.  Figure  4  shows  the  beam 
profiles  at  2  cm  depth  in  the  solid  water  phantom.  It  can  be  seen  that  for  a  20  MeV  electron 
beam,  0.5  cm  leaf  shapes  are  still  distinguishable  on  the  surface  but  become  very  blurred  at 
2  cm  depth.  For  a  6  MeV  electron  beam,  however,  the  effect  of  electron  scattering  becomes  so 
severe  that  a  leaf  width  smaller  than  1.0  cm  will  not  result  in  any  improvement  in  the  spatial 
resolution.  However,  a  small  leaf  width  may  have  the  advantage  of  defining  the  field  more 
precisely  in  the  direction  perpendicular  to  the  leaves. 

Based  on  these  experimental  results,  we  further  performed  Monte  Carlo  simulations 
of  electron  fields  collimated  by  1  cm  wide  leaves  to  study  the  effect  of  material  type  and 
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Figure4.  Beam  profiles  measured  at  2  cm  depth  in  a  solid  water  phantom  for  the  6  MeV  and  20MeV 
electron  fields  shown  in  figure  2:  {a)  20  MeV  along  A-A;  (h)  20  MeV  along  B-B;  (c)  6  MeV  along 
A-A;  and  (d)  6  MeV  along  B-B. 


leaf  thickness.  Although  the  beam  penumbral  widths  did  not  change  significantly  for  leaf 
thicknesses  smaller  than  2  cm,  the  beam  intensity  outside  the  field  was  affected  by  the  leaf 
thickness  and  the  atomic  number  of  the  leaf  material.  As  shown  in  figure  5  for  a  20  MeV  electron 
beam,  1 .5  cm  thick  zinc  reduced  the  electron  fluence  outside  the  field  to  about  5%  of  the  central 
axis  value  (figure  5(a)).  These  electrons  were  mainly  generated  by  the  bremsstrahlung  photons 
in  the  MLC  leaves.  This  was  confirmed  by  the  photon  fluence  as  shown  in  figure  5(^),  where 
1,5  cm  zinc  MLC  leaves  resulted  in  about  60%  higher  photon  fluence  outside  the  field  compared 
with  the  central  axis  photon  fluence.  Some  electrons  were  also  scattered  off  the  leaf  ends  and 
by  air.  For  1 .5  cm  copper,  1.5  cm  lead  and  2  cm  steel,  the  electron  fluence  was  about  2,5%  of 
the  central  axis  value.  The  electron  fluence  was  reduced  to  about  1 .5%  if  the  leaves  were  made 
of  1.5  cm  tungsten.  This  was  reflected  by  the  30%  smaller  photon  fluence  under  the  tungsten 
MLC  leaves  compared  with  the  central  axis  photon  fluence.  Clearly,  tungsten  was  superior  to 
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Figure  5.  Monte  Carlo  simulated  electron  (a)  and  photon  {b)  planar  fluence  in  the  penumbral 
region  and  outside  the  treatment  field  for  a  Varian  Clinac  2100C  20  MeV  electron  beam  collimated 
by  an  electron  MLC  of  different  leaf  materials  and  thicknesses.  The  air  gap  between  the  electron 
MLC  and  the  scoring  plane  is  7  cm. 


Other  materials  in  terms  of  leaf  leakage.  If  we  increased  the  tungsten  leaf  thickness  to  2  cm 
the  electron  fluence  would  be  reduced  to  less  than  1%  of  the  central  axis  value  and  the  photon 
leakage  would  be  reduced  to  about  50%  of  the  central  axis  value  (not  shown). 

To  study  the  overall  effect  of  the  leaf  leakage,  leaf  scattering,  air  scattering  and  the  extended 
source  in  an  electron  beam,  we  compared  the  dose  distributions  for  single  fields  and  multiple 
abutting  fields  collimated  by  an  electron  MLC  with  1.5  cm  thick  tungsten  leaves.  Figure  6 
shows  the  Monte  Carlo  calculated  dose  distributions  for  a  single  4  cm  x  4  cm  electron  field 
and  a  multiple  abutting  field  of  the  same  size  formed  by  four  1  cm  x  4  cm  electron  fields.  For 
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Figure  6.  Monte  Carlo  calculated  dose  distributions  in  a  water  phantom  for  Varian  Clinac  6  and 
20  MeV  electron  beams  collimated  by  an  electron  MLC  of  1 .5  cm  thick  tungsten  leaves  for  a  single 
4  cm  X  4  cm  electron  field  and  a  4  cm  x  4  cm  field  formed  by  four  1  cm  x  4  cm  electron  fields: 
(a)  dose  at  surface  for  a  20  MeV  beam;  (h)  dose  at  3  cm  depth  for  a  20  MeV  beam;  (c)  dose  at 
surface  for  a  6  MeV  beam. 


a  20  MeV  electron  beam,  the  dose  at  the  phantom  surface  for  the  abutting  field  shows  about 
4%  fluctuation  compared  with  a  single  electron  field  (figure  6(a)).  This  is  potentially  due  to 
the  effect  of  leaf  shape  and  extended  source.  The  dose  outside  the  field  for  the  abutting  field  is 
about  three  times  higher  than  that  for  the  single  field,  which  is  mainly  caused  by  the  leaf  leakage 
due  to  the  longer  beam-on  time  to  deliver  the  four  1  cm  x  4  cm  fields  and  electron  scattering 
off  the  leaf  ends.  The  dose  at  3  cm  depth  shows  little  difference  between  the  abutting  field  and 
the  single  field  except  for  the  dose  near  the  field  edges  and  outside  the  field  (figure  6(b)).  For  a 
6  MeV  electron  beam,  the  dose  at  the  phantom  surface  for  the  abutting  field  is  almost  the  same 
as  that  for  the  single  field  (figure  6(c)).  The  dose  outside  the  field  for  the  abutting  field  is  only 
slightly  higher  than  that  for  the  single  field.  The  effect  of  leaf  leakage  is  very  small  for  a  6  MeV 
beam  and  the  dose  immediately  outside  the  field  is  thought  to  be  mainly  due  to  the  effect  of 
electron  scattering  in  the  air.  It  seems  that  field  abutting  with  1  cm  beamlets  collimated  by 
an  electron  MLC  can  provide  adequate  beam  characteristics  for  MERT  for  the  beam  energies 
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investigated.  However,  the  dose  outside  the  field  needs  to  be  minimized  through  beam  energy 
and  leaf  sequence  optimization. 

3.2,  Comparisons  of  a  photon  MLC  and  an  electron  MLC 

There  have  been  studies  on  electron  beam  collimation  using  a  photon  MLC  (Karlsson  et  al 
1999).  One  of  the  advantages  of  using  a  photon  MLC  is  the  possibility  of  combining  both 
photon  and  electron  beams  in  the  same  plan.  An  essential  requirement  for  matching  a  photon 
beam  and  an  electron  beam  at  different  depths  is  that  both  beams  share  the  same  source 
position.  Karlsson  et  al  (1999)  proposed  several  modifications  to  the  design  of  a  Varian 
Clinac  2300CD  accelerator,  one  of  which  was  to  replace  the  intervening  air  with  helium.  This 
could  significantly  reduce  the  effect  of  electron  scattering  in  the  air  on  the  beam  penumbra. 
However,  filling  the  accelerator  head  with  helium  requires  major  modifications  to  the  existing 
accelerator  design.  In  this  work,  we  have  investigated  an  alternative  solution — a  thin  leaf 
MLC  at  the  electron  cutout  level  to  reduce  the  air  scattering  effect.  As  can  be  seen  in  figure  7, 
the  unfocused  MLC  leaf  ends  could  scatter  the  electrons  very  significantly  to  degrade  the 
beam  characteristics  near  the  field  edges.  The  Varian  MLC  has  rounded  leaf  ends,  which 
are  expected  to  have  similar  dosimetric  characteristics  as  the  unfocused  MLC  studied  here. 
Focused  leaf  ends  could  greatly  improve  the  beam  edges  and  provided  even  slightly  better 
dose  profiles  inside  the  field  for  a  20  MeV  electron  beam  compared  with  an  electron  MLC 
(figures  7(a)“(c)),  primarily  due  to  the  reduction  of  electron  scattering  in  the  accelerator  head 
(helium  versus  air).  The  dose  outside  the  field  was  slightly  lower  for  the  electron  MLC  than  for 
the  photon  MLC.  For  a  6  MeV  beam,  an  electron  MLC  gave  slightly  better  surface  dose  profiles 
both  inside  and  outside  the  field  than  the  focused  and  unfocused  photon  MLC.  However,  the 
dose  profiles  became  practically  similar  at  the  depth  of  the  maximum  dose  and  greater  depths 
(not  shown).  Note  that  in  these  comparisons,  we  have  placed  the  phantom  surface  at  20  cm 
below  the  photon  MLC  and  7  cm  below  the  electron  MLC  to  minimize  the  effect  of  electron 
scattering  in  the  air  or  helium  between  the  MLC  and  the  phantom.  It  is  evident  that  an  electron 
MLC  will  have  similar  dosimetric  characteristics  as  a  photon  MLC  with  focused  leaf  ends  but 
without  the  need  to  replace  the  air  in  the  accelerator  head  with  helium. 


3.3.  Comparison  of  beamlet  distributions 

The  accuracy  of  the  beamlet  distribution  calculation  may  play  an  important  role  in  the  treatment 
planning  optimization  process.  Ma  etal(\ 999)  reported  significant  differences  in  the  final  dose 
distributions  of  the  optimized  treatment  plans  computed  by  a  commercial  inverse  treatment 
planning  system  with  a  finite-size  pencil  beam  and  the  Monte  Carlo  method.  Pawlicki  et  al 
(1999)  demonstrated  that  inaccurate  beamlet  distributions  may  result  in  under-dosing  in  the 
target  and  over-dosing  in  the  adjacent  critical  structures,  and  using  the  Monte  Carlo  calculated 
beamlets  could  potentially  reduce  the  uncertainty  in  the  photon  IMRT  dose  distributions.  This 
was  demonstrated  again  by  Jeraj  and  Keall  (1999)  using  a  Monte  Carlo  dose  calculation  based 
inverse  planning  algorithm. 

It  has  been  shown  that  the  electron  beam  dose  distributions  calculated  by  the  pencil  beam 
algorithm  as  implemented  in  some  commercial  treatment  planning  systems  could  be  fairly 
uncertain  in  the  regions  near  material  interfaces  and  inhomogeneities  (Cygler  et  al  1987, 
Shorn  et  al  1986,  Mackie  et  al  1994,  Ma  et  al  1999).  We  have  computed  the  beamlet  dose 
distributions  using  the  3D  pencil  beam  as  implemented  in  the  FOCUS  treatment  planning 
system  (Computerized  Medical  Systems,  St  Louis,  MO)  and  compared  them  with  the  Monte 
Carlo  calculated  beamlets.  Figure  8  shows  the  dose  distributions  calculated  using  the  Monte 
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Figure?.  Monte  Carlo  calculated  dose  distributions  in  a  water  phantom  for  a  10  cm  x  10  cm  field 
collimated  by  an  electron  MLC  with  1 .5  cm  thick  tungsten  leaves  and  a  photon  MLC  with  7.5  cm 
thick  leaves  on  a  Varian  Clinac  2100C  accelerator:  (a)  surface  dose  for  a  20  MeV  beam;  (b)  dose 
at  3  cm  depth  for  a  20  MeV  beam;  (c)  dose  at  6  cm  depth  for  a  20  MeV  beam;  (d)  surface  dose  for 
a  6  MeV  beam.  The  electron  MLC  has  straight  leaf  ends.  The  photon  MLC  has  either  straight  or 
double-focused  leaves. 


Carlo  method  (a,  c,  e)  and  the  FOCUS  3D  pencil  beam  algorithm  (b,  d,  f)  for  a  1  cm  x  I  cm 
12  MeV  beamlet  incident  on  a  patient  phantom  built  from  CT  data.  For  beamlets  with  normal 
incidence  (figures  8(«)  and  (/?)),  the  difference  in  the  dose  distributions  in  the  heart  was  evident: 
the  Monte  Carlo  calculated  isodose  lines  varied  with  the  heart  contour  while  the  pencil  beam 
isodose  lines  remained  symmetrical  despite  the  change  in  material  densities.  Figures  8(c) 
and  {d)  show  the  beamlet  distributions  with  a  10  cm  air  gap.  The  difference  is  clearly  seen 
near  the  surface.  The  beamlet  distributions  again  differed  significantly  in  the  lung  for  oblique 
incidence  (figures  8(c)  and  (/*)).  The  axis  of  the  beamlet  was  intentionally  placed  to  go  through 
soft  tissues  and  bones.  The  pencil  beam  isodose  lines  seemed  to  stretch  according  to  the  beam 
axis  pathlength  and  showed  no  signs  of  electron  build-down  near  the  low-density  material. 
These  results  provided  enough  evidence  to  show  that  to  ensure  the  accuracy  of  the  optimized 
dose  distributions  for  MERT  we  should  use  the  Monte  Carlo  method  to  compute  the  electron 
beamlets  for  the  inverse  planning  process. 


Modulated  electron  beam  for  radiotherapy 


2305 


e  f 

Figures.  Dose  distributions  calculated  using  the  Monte  Carlo  method  (a,  c,  e)  and  the  FOCUS 
3D  pencil  beam  algorithm  {b,  d,  /)  for  a  1  cm  x  1  cm  12  MeV  beamlet  with  normal  incidence 
{a  and  h),  normal  incidence  plus  10  cm  air  gap  (c  and  d),  and  oblique  incidence  (<?  and  /).  The 
beamlet  size  is  defined  at  100  cm  SSD.  The  isodose  lines  shown  are  10,  20,  30, 40, 50, 60, 70,  80 
and  90%  of  the  maximum  dose  respectively. 


3 A.  MERT  versus  photons:  a  hypothetical  treatment  plan 

Modulated  electron  radiotherapy  is  a  general  purpose  technique  that  should  be  advantageous 
in  many  clinical  situations.  An  exhaustive  investigation  of  the  specific  advantages  of  MERT 
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Figure  .  Treatment  plans  for  a  hypothetical  breast  case  using  tangential  6  MV  photon  beams  {a) 
and  MERT  with  6, 12  and  20  MeV  electron  beams  {b).  Both  plans  were  calculated  using  the  Monte 
Carlo  method.  The  isodose  lines  (90, 70,  50  and  30%)  represented  50,  38.9, 25  and  16.7  Gy. 


over  traditional  treatment  modalities  on  a  site  by  site  basis  is  outside  the  scope  of  this  work. 
However,  to  demonstrate  the  possibility  of  improving  the  dose  homogeneity  in  the  target  and  the 
reduction  of  the  dose  in  surrounding  normal  tissues,  we  compare  the  dose  distributions  to  treat 
a  hypothetical  target  using  tangential  photon  beams  and  MERT.  The  purpose  of  the  comparison 
was  to  illustrate  the  concept  of  MERT  but  not  to  draw  specific  conclusions  on  the  use  of  either 
technique.  Previous  investigators  have  used  artificial  phantoms  and  hypothetical  targets  to 
mimic  different  treatment  sites  (e.g.  Hyodynmaa^f  a/  1996,  Asell  etal  1997,  Ebert  and  Hoban 
1997).  We  considered  it  to  be  clinically  relevant  to  use  a  more  realistic  patient  geometry  (built 
from  CT  data)  in  our  comparison,  although  the  target  definition  and  beam  setup  are  somewhat 
arbitrary.  Figure  9  shows  the  hypothetical  treatment  plan  using  tangential  6  MV  photon  beams 
and  MERT  with  normally  incident  6,  12  and  20  MeV  electron  beams.  The  intensity  maps  for 
each  electron  beam  energy  are  shown  in  figure  10.  The  beamlet  size  was  1  cm  x  1  cm  at  1 00  cm 
at  isocentre.  It  is  worth  noting  here,  that  as  a  matter  of  practicality,  it  is  impossible  to  create 
the  intensity  maps  shown  in  figure  10  using  the  conventional  electron  cutout  approach  but  the 
electron  MLC  is  a  viable  alternative.  The  dose  distributions  for  both  plans  were  calculated 
using  the  Monte  Carlo  method.  The  isodose  lines  (90,  70,  50  and  30%)  were  normalized  in 
such  a  way  that  the  90%  isodose  surface  would  receive  the  prescribed  target  dose  of  50  Gy. 
For  the  photon  plan,  the  90%  dose  line  also  included  a  margin  in  the  lung  to  account  for  the 
effect  of  patient  breathing.  This  was  not  needed  for  the  electron  plan  as  the  electron  beams 
were  incident  enface  and  the  electron  beamlet  dose  distributions  do  not  vary  significantly  with 
breathing.  Figure  1 1  shows  the  dose  volume  histograms  (DVH)  for  the  hypothetical  treatment 
plans  shown  in  figure  9.  The  target  DVH  together  with  the  right  lung  DVH  are  shown  in 
figure  11  (fl)  (as  percentage  volume)  and  the  right  lung  DHV  and  the  ‘total  body’  (including 
everything  inside  the  external  contour)  DVH  are  shown  in  figure  1 1(^)  (as  absolute  volume). 

It  is  clear  that  MERT  provided  better  dose  homogeneity  in  the  target  region  than  tangential 
photon  beams.  Tangential  photon  beams  produced  hot  spots  in  the  target  and  cold  spots  near 
the  skin  (figure  1 1(«)).  MERT  significantly  reduced  the  dose  to  the  lung  relative  to  tangential 
photon  beams;  the  maximum  dose  to  the  lung  was  reduced  from  50  Gy  for  a  tangential  treatment 
to  35  Gy  for  MERT  (figure  1 1(a)).  However,  MERT  increased  the  volume  of  the  lung  that 
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Figure  10.  Intensity  maps  for  the  three  electron  beam  energies.  The  beamlet  size  for  each  port 
was  1  X  1  cm“.  Darker  beamlets  indicate  a  higher  weight  than  the  lighter  beamlets  and  the  grey 
scale  for  all  three  maps  is  in  absolute  terms. 


received  a  lower  dose  (10%  more  volume  received  5  Gy  and  20%  more  volume  received  2  Gy) 
compared  with  tangential  photon  beams.  The  clinical  significance  of  the  increased  lung  volume 
receiving  such  a  low  dose  needs  to  be  investigated.  On  the  other  hand,  over  150  cm^  of  lung 
received  much  less  dose  with  MERT  compared  with  tangential  photon  beams,  which  could 
result  in  reduced  lung  complications  (figure  11(^)).  Another  clear  benefit  with  MERT  is  the 
exclusion  of  the  surrounding  normal  tissue  from  the  high  dose  volume  (figure  1 1(^)).  Over 
1000  cm^  of  normal  tissue  received  10-30  Gy  less  dose  in  a  MERT  plan  compared  with  this 
tangential  photon  beam  plan. 

4.  Conclusions 

In  this  work,  we  have  investigated  the  feasibility  of  modulating  both  energy  and  intensity  of 
electron  beams  for  radiotherapy.  This  was  achieved  by  combining  electron  beams  of  different 
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Dose  (Gy) 
ib) 


Figure  11.  Dose  volume  histograms  (DVH)  for  the  breast  treatment  plans  shown  in  figure  9: 
(a)  DVH  shown  as  percentage  volume  for  the  target  (PTV)  and  the  right  lung  and  (fr)  DVH  shown 
as  absolute  volume  for  the  right  lung  and  the  ‘total  body’  which  includes  everything  inside  the 
external  contour. 


nominal  energies  and  variable  intensity  distributions.  A  prototype  electron  MLC  was  built  to 
study  the  characteristics  of  MLC-colIi mated  electron  beams  and  the  Monte  Carlo  simulations 
were  used  to  investigate  the  effect  of  MLC  leaf  material,  thickness,  shape  and  location.  The 
beamlet  distributions  calculated  using  a  3D  electron  pencil  beam  algorithm  as  implemented 
in  a  commercial  treatment  planning  system  and  the  Monte  Carlo  method  were  compared  for 
electron  beams  of  different  energies,  extended  air  gaps,  oblique  incidence  and  heterogeneous 
geometries.  A  hypothetical  breast  case  was  used  to  compare  the  dose  distributions  using 
tangential  photons  and  MERT  for  target  coverage  (dose  homogeneity)  and  normal  tissue  sparing 
(dose  reduction  in  the  lung  and  other  surrounding  normal  tissues). 
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Our  results  showed  that  an  electron  MLC  at  the  electron  cutout  location  can  provide 
adequate  beam  colli mation  for  MERT  without  the  need  to  replace  the  air  in  the  accelerator 
head  and  between  the  MLC  and  the  phantom  with  helium.  The  beam  characteristics  collimated 
by  an  electron  MLC  are  comparable  to  those  collimated  by  a  focused  photon  MLC.  However, 
the  latter  requires  the  accelerator  head  and  between  the  MLC  and  the  phantom  to  be  filled  with 
helium,  which  may  be  impractical  for  some  accelerators  because  of  the  major  modifications 
needed  to  the  structure  design.  An  electron  MLC  can  also  be  used  in  place  of  a  cutout.  The 
Monte  Carlo  method  can  accurately  simulate  particle  transport  in  cases  involving  extended  air 
gaps,  oblique  incidence  and  heterogeneous  anatomy,  and  is  therefore  suitable  for  the  beamlet 
calculation  for  MERT  treatment  optimization.  Our  preliminary  results  based  on  a  hypothetical 
breast  case  demonstrated  the  potential  of  MERT  for  uniform  target  coverage  and  normal  tissue 
sparing.  To  fully  explore  the  potential  of  MERT,  further  studies  need  to  be  carried  out  for 
realistic  clinical  cases  and  for  other  treatment  sites  such  as  the  head  and  neck. 
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Abstract  A  new  EGS4/PRESTA  Monte  Carlo  user  code,  MCDOSE,  has  been  developed  as 
a  routine  dose  calculation  tool  for  radiotherapy  treatment  planning.  It  is  suitable  for  both 
conventional  and  intensity  modulated  radiation  therapy.  Two  important  features  of  MCDOSE 
are  the  inclusion  of  beam  modifiers  in  the  patient  simulation  and  the  implementation  of  several 
variance  reduction  techniques.  Before  this  tool  can  be  used  reliably  for  clinical  dose  calculation,  it 
must  be  properly  validated.  The  validation  for  beam  modifiers  has  been  performed  by  comparing 
the  dose  distributions  calculated  by  MCDOSE  and  the  well-benchmarked  EGS4  user  codes  BEAM 
and  DOSXYZ.  Various  beam  modifiers  were  simulated.  Good  agreement  in  the  dose  distributions 
was  observed.  The  differences  in  electron  cutout  factors  between  the  results  of  MCDOSE  and 
measurements  were  within  2%.  The  accuracy  of  MCDOSE  with  various  variance  reduction 
techniques  was  tested  by  comparing  the  dose  distributions  in  different  inhomogeneous  phantoms 
with  those  calculated  by  DOSXYZ  without  variance  reduction.  The  agreement  was  within  1,0%. 
Our  results  demonstrate  that  MCDOSE  is  accurate  and  efficient  for  routine  dose  calculation  in 
radiotherapy  treatment  planning,  with  or  without  beam  modifiers. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version;  see  www .  iop .  org) 


1.  Introduction 

With  the  rapid  development  of  computer  technology,  the  Monte  Carlo  technique,  currently 
the  most  accurate  method  for  dose  calculation  (Shortt  et  al  1986,  Mackie  1990,  Rogers  and 
Bielajew  1990,  Andreo  1991,  Mackie  et  al  1994,  DeMarco  et  al  1998,  Ma  et  al  1999a),  is 
becoming  more  practical  for  use  in  radiation  therapy  treatment  planning  systems.  It  has  been 
argued  that  the  outcome  of  radiation  therapy  treatment  may  be  improved  by  using  Monte  Carlo 
dose  calculation  (Mohan  1997,  Nahum  1997),  especially  for  intensity  modulated  radiotherapy 
(IMRT)  treatment  planning  (Jeraj  and  Keall  1999,  Ma  etal  1999a,  2000a,  b).  The  Monte  Carlo 
method  is  being  developed  for  treatment  planning  dose  calculations  by  several  groups  (Mackie 
1990,  Hartmann-Siantar  et  al  1997,  Ma  et  al  1997,  1999a,  DeMarco  et  al  1998,  Wang  et  al 
1998,  Mubata  et  al  1998,  Libby  etal  1998,  Faddegon  et  al  1998). 

At  present,  there  are  several  general-purpose  Monte  Carlo  codes  in  widespread  use 
for  radiation  transport  simulation,  e.g.  Electron  Gamma  Shower  version  4  (EGS4)  (Nelson 
et  al  1985),  ETR AN/ITS  (Berger  and  Seltzer  1987)  and  Monte  Carlo  A^-particle  (MCNP) 
(Hendricks  1994)  and  PENELOPE  (Salvat  et  al  1996).  It  is  well  known  that  the  EGS4  code 
system  is  very  well  documented  and  it  has  been  thoroughly  benchmarked  in  the  energy  region 
of  dosimetric  interest  (Rogers  1984a,  b,  Rogers  and  Bielajew  1984,  1986,  1990).  DOSXYZ  is 
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an  EGS4  based  Monte  Carlo  simulation  code  for  calculation  dose  distribution  in  a  rectilinear 
voxel  phantom  (Rogers  etal  1995,  MdiCtal  1995).  It  is  also  well  benchmarked.  Comparisons 
of  the  DOSXYZ  results  with  measurements  have  been  reported  previously  (Rogers  et  al  1995, 
Kapur  et  al  1998,  Ma  1998,  Zhang  et  al  1999,  Ma  et  al  1999a,  2000a).  Unfortunately,  these 
codes  are  too  slow  to  be  acceptable  for  routine  dose  calculation  in  treatment  planning  systems 
with  existing  hardware.  Some  codes  based  on  the  Monte  Carlo  method  have  been  developed 
to  speed  up  the  dose  calculation,  such  as  Macro  Monte  Carlo  (MMC)  (Neuenschwander  and 
Bom  1992,  Neuenschwander  et  al  1995),  Voxel  Monte  Carlo  (VMC  and  XVMC)  (Kawrakow 
et  al  1996,  Fippel  1999)  and  Super  Monte  Carlo  (SMC)  (Keall  and  Hoban  1996a).  Further 
verification  and  improvements  are  needed  for  routine  treatment  planning  dose  calculations. 

A  new  EGS4/PRESTA  (Bielajew  and  Rogers  1987)  user  code,  MCDOSE,  has  been 
developed  for  routine  clinical  dose  calculation  (Ma  et  al  1999b,  2000c).  MCDOSE  was 
designed  as  a  dose  calculation  module  for  easy  implementation  in  a  radiotherapy  treatment 
planning  system.  We  have  implemented  it  in  an  existing  commercial  RTF  system  for 
conventional  photon/electron  beams  (Ma  et  al  1999a)  and  intensity  modulated  radiotherapy 
(IMRT)  dose  verification  (Ma  etal  2000a).  MCDOSE  shares  similar  structure  and  geometry 
definition  with  other  EGS4  user  codes,  such  as  DOSXYZ.  MCDOSE  and  DOSXYZ  have 
some  common  features.  They  can  both  simulate  the  transport  of  photons  and  electrons  (and 
positrons)  in  a  3D  rectilinear  phantom  geometry.  The  volume  elements  (voxels)  in  the  phantom 
can  have  uniform  or  variable  dimensions  and  the  material  in  the  voxel  can  be  specified  by  the 
user  or  determined  from  the  electron  density  data  derived  from  the  patient  CT  data.  Source 
models  are  supported  in  both  codes.  But  there  are  other  features  in  MCDOSE,  which  are 
absent  in  DOSXYZ: 

(a)  Advanced  multiple-source  models  (Ma  et  al  1997,  Ma  1998,  Deng  et  al  2000,  Jiang  et  al 
2000a,  b)  are  used  as  source  input  for  both  photon  and  electron  beams  in  MCDOSE.  The 
multiple-source  models  implemented  in  MCDOSE  are  more  accurate  and  flexible,  and 
easy  for  commissioning. 

(b)  Beam  modifiers  such  as  jaws,  wedges,  blocks,  compensators,  electron  cutouts  and  bolus 
are  included  by  MCDOSE  in  the  patient  simulation.  Both  static  and  dynamic  MLC  fields 
can  be  simulated  for  conventional  and  intensity  modulated  radiotherapy. 

(c)  Several  variance-reduction  techniques  have  been  implemented  in  MCDOSE.  Both  codes 
have  range  rejection.  However,  it  was  implemented  in  different  ways  and  in  different 
simulation  procedures. 

(d)  There  is  an  option  for  selecting  the  geometry  coordinates  from  two  definitions  in 
MCDOSE.  One  (x,y,z)  is  the  same  as  that  in  DOSXYZ  for  phantom  study.  The  other 
{x\  y ,  zO  is  for  convenient  patient  dose  calculation  corresponding  to  the  treatment  system. 
The  relationship  between  the  two  coordinates  is  jc  =  x',  y  =  and  z  =  y'. 

(e)  Beamlets  dose  calculation  for  Monte  Carlo  inverse  planning  can  be  performed  by 
MCDOSE  for  both  photon  and  electron  beams. 

(f)  MCDOSE  can  produce  dose  volume  histograms  (DVH)  using  the  patient  contour 
information. 

The  inclusion  of  beam  modifiers  and  the  implementation  of  variance  reduction  techniques 
makes  MCDOSE  more  practical  for  routine  dose  calculation.  Beam  modifiers  play  an 
important  role  in  radiotherapy  treatment.  They  modify  the  beam  shape  and/or  intensity 
distribution  and  therefore  change  the  dose  distribution  in  the  patient.  They  need  to  be 
considered  for  treatment  planning  dose  calculations.  It  is  not  sufficient  to  treat  jaws  and  blocks 
as  apertures  and  to  change  the  profile  using  the  wedge  or  compensator  factors  conventionally. 
The  effect  of  beam  modifiers  has  been  studied  by  Schach  von  Wittenau  et  al  (2000).  It  was  found 
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that  the  scattered  electrons  generated  in  jaws  could  increase  the  surface  dose  by  several  per  cent 
depending  on  the  beam  energy  and  field  size.  The  photon  transmission  and  bremsstrahlung 
generated  in  the  modifiers  need  to  be  considered.  To  ensure  simulation  accuracy,  these  beam 
modifiers  should  be  simulated  with  accurate  geometry  when  the  Monte  Carlo  method  is  used 
for  routine  dose  calculation.  Some  codes,  like  BEAM  and  PEREGRINE  (Walling  et  al  1998), 
can  be  used  to  simulate  beam  modifiers.  But  BEAM  can  only  simulate  them  separately  from 
patient  dose  calculation.  One  has  to  simulate  particle  transport  in  the  beam  modifiers  using 
BEAM  first.  After  scoring  the  phase  space  below  the  modifiers,  dose  calculation  in  the  patient 
phantom  can  then  be  performed  (Ma  et  al  1999a).  This  two-step  method  works  well  but  is 
inconvenient  and  can  be  time-consuming.  It  is  not  suitable  for  routine  dose  calculations.  The 
inclusion  of  beam  modifiers  into  the  patient  dose  calculation  simplifies  the  simulation  procedure 
in  terms  of  simulation  steps  and  intermediate  data  storage.  Although  Monte  Carlo  is  an  accurate 
method,  the  CPU  time  needed  to  perform  statistically  meaningful  dose  calculation  is  always  a 
problem  that  prevents  it  from  being  applied  in  routine  dose  calculations.  The  variance  reduction 
techniques  implemented  in  MCDOSE  generally  speed  up  the  dose  calculation  by  a  factor  of 
1 0-30  (Maer^/ 2000c).  This  makes  the  code  more  practicable  for  routine  application.  The  dose 
calculation  using  MCDOSE  for  a  nine-field  IMRT  plan  (including  pre-  and  post-optimization 
dose  calculation)  can  be  finished  in  1  to  4  h  on  a  450  MHz  Pentium  III  personal  computer 
(Pawlicki  et  al  1999,  Ma  et  al  2000b). 

Before  MCDOSE  can  be  used  reliably  for  clinical  dose  calculation,  it  must  be  extensively 
validated.  The  accuracy  of  the  dose  calculation  depends  on  the  input  beam  data  and  the 
implementation  of  the  user  code.  Source  modelling  and  beam  commissioning  have  been 
investigated  by  Ma  et  al  (1997),  Ma  (1998),  Deng  et  al  (2000)  and  Jiang  et  al  (2000a,  b). 
Excellent  agreement  in  dose  distributions  in  both  homogeneous  and  heterogeneous  phantoms 
has  been  achieved  between  source  model,  phase  space  data  and  measurements.  In  this  work,  we 
focus  on  the  validation  of  the  user  code  for  beam  modifiers  and  variance  reduction  techniques. 
The  geometry  description  for  beam  modifiers  and  the  implementation  of  variance  reduction 
techniques  may  affect  the  simulation  accuracy  of  the  code. 

The  validation  of  beam  modifiers  was  performed  by  comparing  the  simulation  results  of 
MCDOSE  with  those  of  BEAM/DOSXYZ.  The  dose  distributions  in  specifically  designed 
homogeneous  and  inhomogeneous  phantoms  with  different  source  inputs  calculated  using 
MCDOSE  and  DOSXYZ  are  presented.  Machine  data  for  two  Varian  accelerators,  Clinac 
2100C  and  2300C/D,  were  measured  as  required  by  the  commissioning  procedure. 


2.  Materials  and  methods 

2.7.  Simulation  of  beam  modifiers  in  MCDOSE 

In  this  work,  we  have  implemented  beam  modifier  simulation  in  MCDOSE.  The  combination 
of  the  modifiers,  the  geometry  parameters  and  materials  can  be  input  by  the  user.  Several 
different  modifiers  can  be  used  together.  Each  of  them  occupies  a  slab.  For  photon  beams  our 
multiple-source  model  covers  all  the  fixed  components  above  the  Jaws  in  the  accelerator  head 
(Deng  etal  2000).  What  we  need  to  consider  for  routine  dose  calculations  are  the  jaws,  wedges, 
blocks,  compensators,  MLC  and  bolus.  For  electron  beams  our  multiple-source  model  covers 
all  the  components  down  to  the  lowest  scraper  of  the  applicator  (Jiang  et  al  2000a).  Only  the 
electron  cutout  and  bolus  need  to  be  considered  in  MCDOSE. 

The  X  and  y  jaws  are  modelled  as  two  pairs.  The  inner  surface  of  each  jaw  focuses  to  the 
target.  The  input  parameters  are  x  1 ,  x2,  y  1 ,  y 2  and  the  material  for  the  jaws.  Here  x  1 ,  x2,  y  1 , 
y2  specify  the  locations  of  the  inner  edges  of  the  jaw  bars.  They  correspond  to  the  field  size 
at  100  cm  SSD.  A  schematic  diagram  for  the  jaws  is  shown  in  figure  \(a). 
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Figure  1.  Schematic  diagram  of  the  beam  modifiers  simulated  in  MCDOSE.  (rt)  X  and  Y  jaws. 
ijb)  Photon  beam  block  with  tray,  (c)  Wedge  and  tray,  {d)  Compensator  and  tray. 


To  simulate  the  treatment  block  and  the  tray,  the  location,  thickness  and  material  of  the 
block  and  the  tray  are  required.  The  location  is  defined  as  the  distance  to  the  isocentre  from  the 
bottom  surface  of  the  block  tray.  The  user  also  needs  to  specify  the  opening  by  the  coordinates 
of  the  vertices  that  are  projected  at  100  cm  SSD.  The  planes  defining  the  inner  surfaces  of  the 
opening  are  all  angled  with  respect  to  the  Z-axis  towards  the  target  as  a  single  focus  point.  The 
block  opening  can  be  of  any  shape.  The  user  needs  to  input  the  points  continuously  around 
the  perimeter  of  the  opening.  The  beam’s  eye  view  and  X-Z  view  of  the  block  and  tray  are 
shown  in  figure  \(b). 

Wedges  can  be  simulated  by  MCDOSE  using  the  geometry  information  about  the 
orientation,  location,  material  and  the  two-dimensional  point  coordinates  to  specify  its  shape 
(see  figure  1(c)).  Full  simulation  of  the  particle  transport  is  performed.  This  approach 
is  similar  to  the  wedge  definition  for  conventional  dose  calculation  algorithms,  except  for 
one  important  point.  For  the  Monte  Carlo  approach,  the  user  enters  the  exact  wedge 
specifications  without  any  modifications  or  adjustments  in  the  parameters  so  that  the  final 
dose  calculation  matches  the  measured  dose  profiles  of  the  wedge  under  consideration.  Thus, 
in  our  Monte  Carlo  approach,  the  physical  transport  of  particles  in  the  beam  modifiers  and  the 
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reproduction  of  the  phase  space  must  be  accurate  to  obtain  good  agreement  with  measured 
results. 

The  compensator  including  the  tray  can  be  simulated  by  MCDOSE.  The  location,  material 
and  the  tray  thickness  are  required.  The  compensator  is  divided  into  different  bins  in  two 
directions.  The  coordinates  and  thickness  for  each  bin  (a  pixel)  are  required.  An  X-Z  view 
of  the  compensator  and  tray  is  shown  in  figure  \(d).  This  module  can  also  be  used  as  an  MLC 
after  setting  the  tray  material  to  air. 

An  electron  cutout  is  simulated  in  a  way  similar  to  a  photon  block  except  that  it  does  not 
have  a  tray.  The  option  to  have  a  diverging  or  straight  and  parallel  inner  planes  to  define  the 
opening  is  considered  for  different  clinical  applications.  The  location,  material  and  coordinates 
of  each  vertex  to  specify  the  opening  are  required. 

The  bolus  is  simulated  by  adding  an  extra  layer  of  material  to  the  patient’s  geometry 
(outside  the  patient  external  contour)  according  to  its  material  and  thickness. 

2.2.  Simulation  of  beam  modifiers  by  EGS4/BEAM 

In  order  to  validate  the  code  for  beam  modifier  simulation  in  MCDOSE,  EGS4/BEAM  was  used 
to  simulate  block,  wedge,  compensator  and  electron  cutout.  The  simulations  started  with  the 
phase  space  above  the  jaws.  No  variance  reduction  techniques  were  used  in  these  simulations. 
The  cut-off  energy  of  electron  and  photon  were  set  to  be  700  keV  and  10  keV  (total  energy). 
After  the  phase  space  data  were  obtained  below  the  modifiers,  DOSXYZ  was  used  to  carry 
out  the  dose  calculation.  Comparison  between  them  can  prove  that  our  implementation  for  the 
simulation  of  beam  modifier  is  right. 

Two  component  modules,  BLOCK  and  SLABS,  were  used  to  simulate  the  photon  blocks 
and  the  trays.  To  simulate  the  wedge,  the  left  bars  of  two  or  more  pairs  of  JAWS  were  piled 
up  to  match  the  shape.  The  component  modules  MLC  and  SLABS  were  used  to  simulate  the 
compensators  and  trays.  The  component  BLOCK  was  used  to  simulate  an  electron  cutout  with 
the  phase  space  above  the  last  scraper  as  source  input.  Five  subregions  were  used  to  specify 
a  butterfly  shape  for  the  opening  that  was  simulated  for  an  electron  cutout  in  this  paper.  The 
phase  space  data  were  scored  below  the  modifiers.  Then  dose  calculations  were  performed 
using  DOSXYZ  in  a  water  phantom. 

2.3.  Variance  reduction  techniques  and  dose  calculation  details 

The  goal  of  MCDOSE  is  to  perform  quick  and  accurate  dose  calculation  for  radiation 
therapy  treatment  planning.  Several  variance  reduction  techniques  have  been  implemented 
in  MCDOSE  to  speed  up  the  calculation  without  losing  accuracy.  These  techniques  include 
electron  range  and  region  rejection,  photon  interaction  forcing,  particle  splitting,  Russian 
roulette  (Bielajew  and  Rogers  1988)  and  electron  track  repeating  (Kawrakow  et  al  1996,  Keall 
and  Hoban  1996b).  To  speed  up  the  simulation  in  beam  modifiers,  the  electron  region  rejection 
technique  was  applied  to  the  electrons  that  make  little  contribution  to  the  dose.  For  example, 
when  we  transport  particles  in  jaws  and  blocks,  we  have  an  option  to  transport  electrons  only 
in  a  margin  around  the  opening.  Variable  global  cut-off  energy  (ECUT)  has  been  applied  to 
different  regions.  The  techniques  of  photon  interaction  forcing,  splitting  and  Russian  roulette 
are  well  known  and  are  implemented  in  some  codes,  such  as  BEAM/EGS4  system.  They  were 
well  described  by  Bielajew  and  Rogers  ( 1 988)  and  Rogers  etal  (1995).  Electron  track  repeating 
is  similar  to  a  technique  called  correlated  sampling  (Ma  and  Nahum  1993),  in  which  particle 
tracks  were  repeated  in  different  locations  or  geometries  to  improve  computation  efficiency.  It 
was  described  in  great  detail  by  Kawrakow  etal  (1996)  and  Keall  and  Hoban  (1996b).  The  main 
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Figure  2.  Comparisons  of  the  lateral  dose  profiles  at  different  depths  in  a  water  phantom  between 
MCDOSE  and  BEAM/DOS XYZ.  A  45°  wedge  were  simulated  for  a  15  MV  photon  beam,  the  field 
size  was  10  cm  x  10  cm  defined  at  100  cm  SSD. 


idea  is  to  simulate  an  electron  in  a  standard  phantom  such  as  in  water  or  tissue,  recording  all 
the  information  of  the  electron  transport  tracks  including  the  step  length,  direction,  energy 
deposition,  type  of  collision  and  generated  particles  for  each  step,  then  repeat  the  tracks  for 
other  electrons  with  the  same  energy  starting  from  a  different  location  in  the  patient.  The 
initial  direction  for  each  electron  may  be  different.  The  length  of  every  step  in  the  electron 
track  is  inversely  proportional  to  the  material  density  at  that  location.  Stopping  power  ratio  and 
scattering  power  ratio  for  different  materials  need  to  be  considered  to  scale  the  path  lengths 
and  multiple  scattering  angles.  This  technique  saves  the  sampling  time  for  multiple  scattering 
and  collision  for  all  the  repeated  electrons  but  not  for  the  initial  one.  How  many  times  to  repeat 
the  tracks  depends  on  the  problem.  The  number  may  affect  the  final  efficiency.  Our  goal 
was  to  integrate  this  technique  efficiently  with  other  techniques  mentioned  above  to  produce 
a  fast  Monte  Carlo  code.  For  example,  when  photons  go  into  the  phantom,  they  are  split  and 
forced  to  have  same  interactions  at  different  locations.  Many  electrons  are  generated  and  their 
tracks  are  repeated  at  these  locations.  So  time  is  saved  because  we  do  not  perform  initial 
simulation  for  all  the  electrons.  Russian  roulette  is  applied  to  reduce  the  number  of  scattering 
photons. 

To  ensure  the  accuracy  of  the  dose  calculation  with  MCDOSE,  specially  designed 
inhomogeneous  phantoms,  such  as  layered-lung  or  layered-bone  phantoms,  are  used.  The 
dose  distributions  (depth  dose  curves  and  lateral  dose  profiles)  calculated  by  MCDOSE  and 
DOS  XYZ  for  different  photon  beams  and  electron  beams  are  compared.  DOS  XYZ  was 
selected  for  comparison  with  MCDOSE  because  it  has  already  been  tested  thoroughly  and 
it  is  also  an  EGS4  user  code.  Agreement  between  MCDOSE  and  DOSXYZ  means  that  the 
beam  modifier  implementation  in  MCDOSE  is  correct  and  the  application  of  variance  reduction 
techniques  does  not  affect  the  accuracy.  Since  DOSXYZ  has  been  well  benchmarked  by  other 
investigators  and  our  previous  work  (Rogers  et  al  1995,  Kapur  et  al  1998,  Ma  1998,  Zhang 
et  al  1999,  Ma  et  al  1999a,  2000a),  it  is  not  necessary  to  compare  all  the  MCDOSE  results 
with  measurements.  In  this  work,  both  monoenergetic  photon  beams  and  realistic  clinical 
beams  of  Varian  2300C/D  are  studied.  The  goal  for  use  of  monoenergetic  photon  beams  is  to 
make  sure  that  the  code  works  well  for  monoenergetic  beams  (therefore  it  should  also  work 
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Figure  3.  Dose  distributions  in  a  water  phantom  calculated  by  MCDOSE  (thin  curve)  and 
BE  AM/DOSXYZ  (thick  curve)  for  a  10  cm  x  10  cm  field  of  a  15  MV  photon  beam.  The  beam  was 
modified  by  a  hexagonal  block  and  a  60"’  wedge,  (a)  The  isodose  distribution  in  the  X~Z  plane 
through  the  central  axis  (beams  come  from  the  bottom),  (b)  The  isodose  distribution  in  the  X-Y 
plane  at  a  depth  of  3.5  cm. 


for  polyenergetic  beams).  This  test  was  done  thoroughly  in  the  energy  range  of  interest  for 
dosimetry  with  an  interval  of  1  MeV. 

3.  Results 

3.J.  Validation  of  beam  modifiers 

3.1.1.  A  45°  wedge.  Figure  2  shows  the  lateral  dose  profiles  at  different  depths  in  a  water 
phantom  calculated  using  MCDOSE  and  BEAM/DOSXYZ  for  a  45"^  wedge.  Lead  was  used 
as  the  wedge  material.  The  field  size  was  set  to  10  cm  x  10  cm  at  100  cm  SSD.  The  voxel 
size  in  the  water  phantom  was  1.0  cm  x  0.3  cm  x  0.5  cm.  A  15  MV  photon  beam  from  a 
Varian  Clinac  2300C/D  machine  was  used  as  source  input  with  phase  space  scored  above  the 
jaws.  The  comparison  calculations  were  performed  under  the  same  conditions.  The  agreement 
between  the  results  of  MCDOSE  and  BEAM/DOSXYZ  is  within  1%  of  the  maximum  dose. 
The  CPU  time  comparison  between  MCDOSE  and  BEAM/DOSXYZ  is  shown  in  table  1. 
MCDOSE  is  about  16  times  faster  than  BEAM/DOSXYZ  for  this  simulation. 

3.1.2.  A  hexagonal  block  and  a  6(f  wedge.  Figure  3  shows  the  dose  distributions  calculated  by 
MCDOSE  and  BEAM/DOSXYZ  in  a  water  phantom  with  a  hexagonal  block  and  a  60"^  wedge 
for  a  10  cm  X  10  cm  field  of  a  15  MV  photon  beam.  The  block  was  made  of  Cerrobend  and  its 
thickness  was  7.7  cm.  The  block  tray  was  made  of  PMMA  and  its  thickness  was  0.7  cm.  The 
wedge  was  made  of  lead.  The  voxel  size  in  the  water  phantom  was  0.5  cm  x  0.5  cm  x  0.5  cm. 
Figure  ?>{a)  shows  the  isodose  distribution  in  the  X-Z  plane  through  the  central  axis  and 
figure  3(6)  shows  the  isodose  distribution  in  the  X~Y  plane  at  a  depth  of  3.5  cm.  The  difference 
between  the  results  of  MCDOSE  and  BEAM/DOSXYZ  is  within  1%  of  the  maximum  dose. 
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Figure  4.  (a)  The  isodose  distributions  in  XY  view  at  a  depth  of  2.5  cm  in  a  water  phantom  for 
a  modified  10  cm  x  10  cm  15  MV  photon  beam  using  a  specially  designed  compensator.  The 
thick  curve  is  the  simulation  result  of  BEAM/DOSX YZ.  The  thin  curve  is  the  result  of  MCDOSE. 
{b)  The  lateral  dose  profile  along  the  dashed  line  in  (a).  The  curve  with  symbols  in  (b)  is  the  result 
of  MCDOSE. 

Table  1.  The  CPU  time  (hours)  required  by  MCDOSE  and  BEAM/DOSXYZ  for  the  beam  modifier 
simulation  investigated  in  this  work  on  a  Pentium  III  450  MHz  personal  computer.  The  statistical 
uncertainties  of  the  maximum  doses  are  1%. 


Modifier 

45°  wedge 

60°  wedge/block 

E-shaped 

compensator 

Butterfly-shaped 
electron  cutout 

BEAM  +  DOSXYZ 

MCDOSE 

Ratio 

31.5  +  6.1 

2.3 

16.3 

247.5  +  7.9 

13.0 

19.4 

10.3  +  10.9 

5.3 

4.0 

37.9+1.4 

2.5 

15.7 

The  CPU  time  needed  by  MCDOSE  and  BEAM/DOSXYZ  is  given  in  table  1.  MCDOSE  is 
about  19  times  faster  than  BEAM/DOSXYZ  for  this  simulation. 


Validation  of  a  Monte  Carlo  dose  calculation  tool 


2977 


Figure  5.  (a)  is  the  isodose  distributions  in  XY  view  at  a  depth  of  2,25  cm  in  a  water  phantom  for 
a  modified  10  cm  x  10  cm  12  MeV  electron  beam  using  an  electron  cutout  with  a  butterfly-shaped 
opening.  The  thick  curve  is  the  simulation  result  of  BEAM/DOSXYZ.  The  thin  curve  is  the  result 
of  MCDOSE.  {b)  The  lateral  dose  profile  along  the  dashed  line  in  (a).  The  curve  with  symbols  in 
(h)  is  the  result  of  MCDOSE. 


3.L3,  Compensator  Figure  4  shows  the  E-shaped  compensator  isodose  distributions  {a)  and 
lateral  dose  profiles  {b)  at  a  depth  of  2.5  cm  in  the  water  phantom.  They  were  calculated  using 
MCDOSE  and  BEAM/DOSXYZ  respectively.  The  compensator  was  made  of  copper.  The 
thickness  in  the  area  of  the  compensator  forming  the  ‘E’  shape  was  0.5  cm  and  the  thickness  in 
the  remaining  area  was  2  cm.  The  voxel  size  in  the  water  phantom  was  0.4  cm  x  0.4  cm  x  1 .0  cm. 
Agreement  (within  1%)  can  be  observed  from  the  comparison.  In  the  penumbra  regions,  the 
separation  between  the  results  of  MCDOSE  and  BEAM/DOSXYZ  is  within  1  mm.  The  CPU 
time  needed  by  MCDOSE  and  BEAM/DOSXYZ  is  also  shown  in  table  1,  MCDOSE  is  about 
four  times  faster  than  BEAM/DOSXYZ  for  this  simulation. 

3.L4.  Electron  cutout.  Figure  5  shows  the  isodose  distributions  {a)  and  lateral  dose 
profiles  (/?)  in  a  water  phantom  calculated  using  BEAM/DOSXYZ  and  MCDOSE,  fora  12  MeV 
realistic  electron  beamfroma  Varian  2100C  machine  with  a  butterfly-shaped  cutout.  The  cutout 
was  made  of  Cerrobend  and  its  thickness  was  1.4  cm.  The  voxel  size  in  the  water  phantom 
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Figure  6.  Dose  distribution  in  a  tissue-lung-tissue  phantom  calculated  by  DOSX  YZ  and  MCDOSE 
for  6  MeV,  15  MV  and  2  MeV  photon  beams  with  field  size  of  6  cm  x  6  cm  defined  at  1(X)  cm 
SSD.  The  material  in  the  volume  of  —5.0  cm  to  5.0  cm  along  the  X  direction,  -2.0  cm  to  2.0  cm 
along  the  Y  direction  and  7.0  cm  to  12.5  cm  along  the  Z  direction  is  lung.  («)  The  depth  dose 
curves  along  the  central  axis  for  the  beams,  {h)  The  lateral  dose  profiles  along  the  X -axis  at  depths 
of  3.3  cm,  9.0  cm  and  15.0  cm  for  the  15  MV  photon  beam. 


was  0.5  cm  X  0.5  cm  x  0.3  cm.  The  phase  space  of  an  electron  beam  above  the  last  scraper 
of  a  10  cm  X  10  cm  applicator  was  used  as  the  source  input.  Good  agreement  (within  1%) 
has  been  achieved  between  the  results  of  MCDOSE  and  BEAM/DOS XYZ.  In  the  penumbra 
regions,  the  separation  between  MCDOSE  and  BEAM/DOSXYZ  is  within  1  mm.  The  CPU 
time  needed  by  MCDOSE  and  BEAM/DOSXYZ  is  shown  in  table  1.  MCDOSE  is  about  16 
times  faster  than  BEAM/DOSXYZ  for  this  simulation. 

The  electron  beam  cutout  factors  calculated  using  MCDOSE  are  compared  with  the 
published  data  in  table  2.  The  published  data  were  calculated  using  BEAM/DOSXYZ  with  full 
phase  space  input  and  measured  by  Kapur  et  al  (1998).  The  difference  between  the  MCDOSE 
data  and  those  calculated  with  full  phase  space  by  DOSXYZ  is  within  1.0%.  The  differences 
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Figure  7.  Dose  distributions  in  a  tissue-bone-tissue  phantom  calculated  by  DOSXYZ  and 
MCDOSE  for  6  MeV,  15  MV  and  2  MeV  photon  beams  with  field  size  of  6  cm  x  6  cm  defined  at 
1(K)  cm  SSD.  The  material  in  the  volume  of  —5.0  cm  to  5.0  cm  along  the  X  direction,  —2.0  cm  to 
2.0  cm  along  the  Y  direction  and  7.0  cm  to  10.0  cm  along  the  Z  direction  is  bone,  (a)  The  depth 
dose  curve  along  the  central  axis  for  the  beams,  (b)  The  lateral  dose  profiles  along  the  Y  axis  at 
depths  of  3.3  cm,  8.0  cm  and  1 1.0  cm  for  the  15  MV  photon  beam. 

in  the  cutout  factors  between  MCDOSE  data  and  measurement  is  less  than  2%  except  for  the 
cases  of  20  MeV  electron  beams  with  3  cm  x  3  cm  and  4  cm  x  4  cm  inserts. 


3.2.  Validation  of  dose  distribution  in  heterogeneous  phantoms 

To  test  the  accuracy  of  the  dose  distribution  calculated  using  MCDOSE  with  the  variance 
reduction  techniques,  comparisons  are  made  with  the  results  of  DOSXYZ.  Each  example  has 
been  calculated  with  a  source-surface  distance  of  100  cm  and  field  size  of  6  cm  x  6  cm. 
The  global  electron  cut-off  energy  (ECUT)  was  set  to  700  keV  and  photon  cut-off  energy 
(PCUT)  was  set  to  10  keV  in  MCDOSE  and  DOSXYZ,  The  maximum  fractional  energy  loss 
per  electron  step  (ESTEPE)  was  limited  to  4%  for  both  codes.  The  material  compositions  and 
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Figure  8.  Dose  distributions  in  a  tissue-lung-tissue  phantom  calculated  by  DOSXYZ  and 
MCDOSE  for  a  15  MeV  electron  beam  with  field  size  6  cm  x  6  cm.  The  material  in  the  volume 
of  -5.0  cm  to  5.0  cm  along  the  X  direction,  -2.0  cm  to  2.0  cm  along  Y  direction  and  4.2  cm  to 
6.0  cm  along  the  Z  direction  is  lung.  («)  The  depth  dose  curve  along  the  central  axis,  (b)  The 
lateral  dose  profiles  along  the  T-axis  at  depths  of  3.0  cm,  5.0  cm  and  7.0  cm. 


densities  are  taken  from  ICRU  (1992),  Monoenergetic  point  sources  and  phase  space  files  for 
a  15  MV  photon  beam  of  a  Varian  Clinac  2300C/D  linear  accelerator  were  used.  Since  dose 
per  incident  fluence  is  calculated,  absolute  dose  comparisons  are  made  between  MCDOSE  and 
DOSXYZ. 


3.2,1.  Photon  beams.  Figure  6{a)  shows  the  depth  dose  curves  along  the  central  axis  for 
6  MeV  and  2  MeV  monoenergetic  photon  beams  calculated  by  DOSXYZ  and  MCDOSE  in  a 
tissue  phantom  with  3D  lung  inhomogeneity.  At  depths  from  7  cm  to  12.5  cm,  the  phantom 
contains  a  10  cm  x  4  cm  slab  of  lung  in  the  centre.  The  depth  dose  curves  for  a  realistic 
15  MV  photon  beam  from  a  Varian  Clinac  2300C/D  accelerator  calculated  by  DOSXYZ  and 
MCDOSE  in  the  same  phantom  are  also  shown  in  figure  6{a)  and  the  lateral  dose  profiles  at 
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Figure  9.  Dose  distribution  in  tissue-bone-tissue  phantoms  calculated  by  DOSX  YZ  and  MCDOSE 
for  15  MeV  and  6  MeV  electron  beams  with  field  size  6  cm  x  6  cm.  For  the  15  MeV  electron 
beam,  the  material  in  the  volume  of  -5.0  cm  to  5.0  cm  along  the  X  direction,  —2.0  cm  to  2.0  cm 
along  the  Y  direction  and  2.7  cm  to  4.5  cm  along  the  Z  direction  is  bone.  For  the  6  MeV  electron 
beam,  the  material  in  the  volume  of  —5.0  cm  to  5.0  cm  along  the  X  direction,  —2.0  cm  to  2.0  cm 
along  the  Y  direction  and  0.6  cm  to  1 .6  cm  along  the  Z  direction  is  bone,  (a)  The  depth  dose  curves 
along  the  central  axis  for  two  beams,  {b)  The  lateral  dose  profiles  along  the  Y  axis  at  depths  of 
3.4  cm  and  5.0  cm  for  the  15  MeV  electron  beam. 

different  depths  are  shown  in  figure  6{b).  The  phase  space  for  a  15  MV  photon  beam  was 
generated  by  BEAM  and  contains  over  5  million  particles.  The  \a  statistical  uncertainty  on 
the  dose  values  for  these  curves  is  less  than  0.5%. 

Dose  calculations  were  also  performed  for  photon  beams  in  an  inhomogeneous  phantom 
containing  tissue  and  bone.  A  10  cm  x  4  cm  slab  of  bone  was  placed  between  7  cm  and  10  cm 
depth  in  the  centre  of  the  tissue  phantom.  Figure  7(a)  shows  the  depth  dose  curves  along  the 
central  axis  for  different  photon  beams.  Figure  l{b)  shows  the  lateral  dose  profiles  at  different 
depths  for  a  15  MV  photon  beam. 

Figures  6  and  7  show  good  agreement  in  dose  distributions  between  MCDOSE  and 
DOSXYZ  for  different  photon  beams.  The  absolute  doses  calculated  by  these  two  codes 
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Table  2.  The  electron  cutout  factors  for  various  square  inserts  in  10  x  10  cm^  applicator  for  6,  12 
and  20  MeV  beams  calculated  by  MCDOSE  and  compared  with  the  full  phase  space  calculation 
by  DOSXYZ  and  measurement  of  Kapur  et  al  (1998).  The  values  in  the  parentheses  indicate  the 
per  cent  difference  between  the  MCDOSE  calculated  data  and  those  measured  or  calculated  with 
full  phase  space  by  DOSXYZ. 


Energy  (MeV) 

Insert  (cm  x  cm) 

MCDOSE 

Full  phase  space 

Measurement 

6 

3x3 

0.922 

0.923  (-0.1%) 

0.927  (-0.5%) 

4x4 

0.983 

0.982  (0.1%) 

0.988  (-0.5%) 

8x8 

1.009 

1.005  (0.4%) 

1.003  (0.6%) 

12 

3x3 

0.930 

0.930  (0.0%) 

0.928  (0.2%) 

4x4 

0.965 

0.956  (0.9%) 

0.963  (0.2%) 

8x8 

0.999 

1.002  (-0.3%) 

0.991  (0.8%) 

20 

3x3 

0.964 

0.968  (-0.4%) 

0.993  (-2.9%) 

4x4 

0.989 

0.993  (-0.4%) 

1.011  (-2.2%) 

8x8 

0.999 

0.993  (0.6%) 

1.004  (-0.5%) 

in  tissue,  lung,  bone,  even  at  the  interface  between  tissue  and  bone  and  tissue  and  lung  agree  to 
within  1  %  of  the  maximum  dose.  Further  comparisons  of  dose  calculations  for  monoenergetic 
photon  beams  from  1  MeV  to  20  MeV,  and  also  for  6  MV  and  10  MV  photon  beams  were 
performed.  Similar  agreement  between  MCDOSE  and  DOSXYZ  was  achieved.  In  most  cases, 
the  differences  were  within  the  la  statistical  uncertainty  of  1%. 


3.2.2,  Electron  beams.  We  also  performed  dose  calculations  for  monoenergetic  point  source 
electron  beams  with  field  size  6  cm  x  6  cm.  The  depth  dose  curves  of  a  15  MeV  electron 
beam  calculated  by  MCDOSE  and  DOSXYZ  in  a  tissue-lung-tissue  phantom  are  shown  in 
figure  8(«).  The  lung  material  was  at  a  depth  from  4  to  6  cm.  Figure  %{b)  shows  the  lateral  dose 
profiles  at  different  depths  in  the  phantom.  At  a  depth  from  4.2  to  6.0  cm,  the  lung  material 
was  from  —5  cm  to  5  cm  in  the  X  direction  and  from  —2  cm  to  2  cm  in  the  Y  direction.  The 
remaining  volume  was  tissue.  The  difference  between  the  simulation  results  of  MCDOSE  and 
DOSXYZ  is  within  1%  of  the  maximum  dose. 

The  depth  dose  curves  and  lateral  dose  profiles  at  different  depths  for  a  15  MeV  electron 
point  source  beam  in  a  tissue-bone-tissue  phantom  are  shown  in  figures  9(a)  and  (b).  At  a 
depth  from  2.7  to  4.5  cm,  there  was  a  10  cm  x  4  cm  slab  of  bone  in  the  centre  surrounded  by 
tissue.  The  depth  dose  curve  for  a  6  MeV  electron  beam  in  a  tissue-bone-tissue  phantom  is 
also  shown  in  figure  9(a).  The  bone  material  was  at  a  depth  from  0.6  to  1 .6  cm.  The  difference 
between  the  simulation  results  of  MCDOSE  and  DOSXYZ  is  within  1  %  of  the  maximum  dose. 

3.2.3.  CPU  time  comparison.  The  CPU  times  required  by  MCDOSE  and  DOSXYZ  for  dose 
calculations  shown  in  figures  6  and  9  are  listed  in  table  3  fordifferent  photon  and  electron  beams. 
The  calculations  were  performed  on  a  Pentium  III  450  MHz  PC.  The  statistical  uncertainty 
of  the  maximum  dose  is  1%.  The  ratios  in  the  table  are  the  speed-up  factors  of  MCDOSE 
compared  with  DOSXYZ  for  each  beam.  In  general  MCDOSE  is  20-50  times  faster  than 
DOSXYZ  for  photon  beams  and  5-10  times  faster  than  DOSXYZ  for  electron  beams. 


4.  Summary 

A  new  Monte  Carlo  dose  calculation  tool  has  been  developed  and  validated  for  clinical  dose 
calculation.  The  main  features  of  MCDOSE  include  electron  and  photon  beam  reconstruction 
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Table  3.  The  CPU  time  (hours)  required  by  MCDOSE  and  DOSXYZ  for  the  dose  calculations 
shown  in  figures  6  and  9  on  a  Pentium  III  450  MHz  personal  computer.  The  statistical  uncertainties 
of  the  maximum  doses  are  1%. 


Beam 

2  MeV  photon 

6  MeV  photon 

15  MV  photon 

6  MeV  electron 

15  MeV  electron 

DOSXYZ 

9.44 

2.95 

4.86 

0.78 

0.94 

MCDOSE 

0.25 

0.074 

0.098 

0.08 

0.17 

Ratio 

37.7 

39.8 

49.6 

9.8 

5.6 

by  source  models,  simulation  of  beam  modifiers  together  with  the  patient  geometry  for  dose 
calculation  and  application  of  variance  reduction  techniques. 

Comparisons  of  dose  distribution  with  wedges,  blocks,  compensator  and  electron  cutout 
between  MCDOSE  and  BEAM/DOSXYZ  demonstrated  that  dose  calculation  with  modifiers 
in  MCDOSE  is  accurate,  convenient  and  efficient.  To  compare  MCDOSE  and  DOSXYZ  dose 
distributions,  3D  heterogeneous  phantoms  containing  lung  and  bone  were  used.  Excellent 
agreement  between  the  MCDOSE  and  DOSXYZ  results  has  been  obtained  for  both  photon 
and  electron  beams  of  different  energies.  Comparing  with  DOSXYZ,  the  MCDOSE  code  runs 
generally  20-50  times  faster  for  the  photon  beams  and  5-10  times  faster  for  electron  beams 
investigated  in  this  work  without  beam  modifiers. 
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A  hybrid  approach  for  commissioning  electron  beam  Monte  Carlo  treatment  planning  systems  has 
been  studied.  The  approach  is  based  on  the  assumption  that  accelerators  of  the  same  type  have  very 
similar  electron  beam  characteristics  and  the  major  difference  comes  from  the  on-site  tuning  of  the 
electron  incident  energy  at  the  exit  window.  For  one  type  of  accelerator,  a  reference  machine  can  be 
selected  and  simulated  with  the  Monte  Carlo  method.  A  multiple  source  model  can  be  built  on  the 
full  Monte  Carlo  simulation  of  the  reference  beam.  When  commissioning  electron  beams  from 
other  accelerators  of  the  same  type,  the  energy  spectra  in  the  source  model  are  tuned  to  match  the 
measured  dose  distributions.  A  Varian  Clinac  2100C  accelerator  was  chosen  as  the  reference 
machine  and  a  four-source  beam  model  was  established  based  on  the  Monte  Carlo  simulations.  This 
simplified  beam  model  can  be  used  to  generate  Monte  Carlo  dose  distributions  accurately  (within 
2%/2  mm  compared  to  those  calculated  with  full  phase  space  data)  for  electron  beams  from  the 
reference  machine  with  various  nominal  energies,  applicator  sizes,  and  SSDs.  Three  electron  beams 
were  commissioned  by  adjusting  the  energy  spectra  in  the  source  model.  The  dose  distributions 
calculated  with  the  adjusted  source  model  were  compared  with  the  dose  distributions  calculated 
using  the  phase  space  data  for  these  beams.  The  agreement  is  within  1%  in  most  of  cases  and  2% 
in  all  situations.  This  preliminary  study  has  shown  the  capability  of  the  commissioning  approach  for 
handling  large  variation  in  the  electron  incident  energy.  The  possibility  of  making  the  approach 
more  versatile  is  also  discussed.  ©  2000  American  Association  of  Physicists  in 
Medicine.[S0094-2405i00)03m-5] 

Key  words:  electron  beam,  treatment  planning,  Monte  Carlo  simulation,  beam  commissioning, 
source  modeling 


I.  INTRODUCTION 

Electron  beam  radiation  therapy  is  used  extensively  to  treat 
head  and  neck  cancers  to  avoid  the  irradiation  of  the  spinal 
cord,  and  to  treat  chest  walls  to  limit  the  irradiated  volume  of 
lung.  The  currently  available  commercial  systems  for  elec¬ 
tron  treatment  planning  mostly  utilize  the  Hogstrom  algo¬ 
rithm  as  the  dose  calculation  engine,^  which  is  based  on 
Fermi-Eyges  theory.^’^  Due  to  the  inappropriate  treatment  of 
electron  transport  in  inhomogeneous  phantoms,  large  dis¬ 
crepancies  (10%  or  more)  in  the  dose  distributions  have  been 
observed  between  the  current  analytical  algorithms  and  mea¬ 
surements  or  Monte  Carlo  simulations  in  some  clinical  situ¬ 
ations  where  the  treatment  volumes  encompass  air  cavities 
and  bone."*”^  Accurate  dose  calculation  is  an  important  factor 
for  the  widespread  clinical  use  of  electron  therapy  and  the 
development  of  new  electron  therapy  techniques,  such  as 
electron  beam  or  mixed  electron/photon  beam  intensity 
modulated  therapy,  which  are  expected  to  improve  the  con¬ 
formality  of  the  delivered  dose  distribution  to  the  target  vol¬ 
ume  for  some  disease  sites.^"^^ 

The  Monte  Carlo  method  is  generally  considered  to  be  the 
most  accurate  approach  for  electron  dose  calculation  under 
all  circumstances.^ In  particular,  Monte  Carlo  simulation 
can  handle  electron  multiple  scattering  in  the  presence  of 
inhomogeneities  (such  as  bone  and  air  cavity)  much  more 


accurately  than  any  existing  analytical  dose  models.  The  ne¬ 
cessity  of  accurate  electron  dose  calculation  has  motivated 
many  efforts  to  develop  Monte  Carlo  electron  beam  treat¬ 
ment  planning  systems.^’^’*^"^^  Due  to  the  rapid  development 
of  computer  technology  and  the  employment  of  innovative 
variance  reduction  techniques,  it  is  expected  that  treatment 
planning  systems  utilizing  a  Monte  Carlo  dose  engine  will 
begin  to  serve  in  routine  clinical  practice  in  the  next  few 
years.^’^’^®“^^ 

The  commissioning  procedure  for  a  Monte  Carlo  treat¬ 
ment  planning  system  can  be  different  from  that  for  a  con¬ 
ventional  planning  system,  since  it  requires  more  detailed 
and  accurate  clinical  beam  data.^^  For  example,  the  phase 
space  information  (position,  direction  and  energy)  is  needed 
to  represent  particles  coming  out  of  the  accelerator  treatment 
head.  This  information  is  extremely  difficult,  if  not  impos¬ 
sible,  to  acquire  experimentally,  mainly  due  to  the  very  high 
intensity  of  the  therapeutic  electron  beam.^®  Some  research¬ 
ers  tried  to  extract  the  phase  space  information  from  the 
limited  set  of  measured  dose  data  (such  as  depth  dose  curves 
and  dose  profiles)  by  using  a  simple  beam  model. Although 
the  approach  may  have  great  potential,  at  least  currently  it 
uses  many  approximations  and  the  accuracy  of  the  recon¬ 
structed  phase  space  cannot  be  guaranteed.  The  only  method 
to  obtain  the  accurate  electron  beam  phase  space 
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information  is  to  simulate  the  accelerator  treatment  head  us¬ 
ing  the  Monte  Carlo  method.^^"^"^  An  EGS4  Monte  Carlo 
user  code,  omega  beam,  was  developed  specifically  for  this 
purpose. Using  the  BEAM  code,  the  accelerator  treatment 
head  and  electron  applicator  can  be  simulated  to  yield  a  data 
file  containing  the  phase  space  information  for  tens  of  mil¬ 
lions  of  particles  exiting  the  treatment  head.  The  phase  space 
data  can  then  be  used  as  input  to  calculate  dose  distributions 
in  a  patient’s  CT  phantom.^"^ 

However,  direct  simulation  of  the  accelerator  treatment 
head  using  the  Monte  Carlo  method  is  not  a  viable  commis¬ 
sioning  approach  for  Monte  Carlo  treatment  planning.  The 
beam  characteristics  are  usually  different  due  to  variation  in 
accelerator  designs  and  on-site  beam  tuning.  The  simulated 
electron  beam  phase  space  for  one  accelerator  may  not  be 
used  directly  for  another.  It  is  necessary  to  simulate  each 
accelerator  individually  to  obtain  the  phase  space  informa¬ 
tion.  This  fact  presents  three  problems  for  the  clinical  accep¬ 
tance  of  Monte  Carlo  treatment  planning  systems.  First,  the 
simulation  of  the  accelerator  treatment  head  for  every 
energy/applicator  combination  takes  much  more  time  than 
the  commissioning  of  a  conventional  electron  planning  sys¬ 
tem.  As  estimated  by  Faddegon  et  even  for  users  with 
Monte  Carlo  simulation  experience,  it  takes  about  two 
months  of  CPU  time  to  generate  a  complete  set  of  beam  data 
for  a  single  accelerator.  Second,  the  storage  of  the  phase 
space  information  requires  a  lot  of  computer  disk  space.  For 
each  energy/applicator  combination,  a  phase  space  file  is 
usually  pre-calculated  and  stored  in  the  treatment  planning 
computer.  For  accurate  treatment  planning,  a  phase  space  file 
occupies  hundreds  of  megabytes  of  disk  space.  This  is  cer¬ 
tainly  a  substantial  burden  for  the  computer  resources  at 
most  clinical  centers.  Third,  the  generation  and  quality  assur¬ 
ance  of  the  phase  space  data  files  by  simulating  the  treatment 
head  requires  Monte  Carlo  simulation  experience.  Therefore, 
it  is  a  prohibitive  task  for  general  users  to  perform  Monte 
Carlo  simulations  for  their  own  accelerators. 

In  this  paper,  a  hybrid  approach  for  commissioning  elec¬ 
tron  beam  Monte  Carlo  treatment  planning  systems  is  pro¬ 
posed.  This  method  combines  the  advantages  of  the  full 
Monte  Carlo  simulation  and  the  method  of  Janssen  et  It 
is  based  on  the  assumption  that  accelerators  of  the  same  type 
have  very  similar  electron  beam  characteristics  and  the  major 
difference  is  the  electron  incident  energy  at  the  exit  window 
due  to  beam  tuning  during  linac  acceptance.  By  simulating  a 
reference  accelerator  for  a  particular  type  of  accelerator  us¬ 
ing  the  Monte  Carlo  BEAM  code,^"^  a  beam  model  is  con¬ 
structed  using  the  resultant  phase  space  information.  The 
beam  model  is  a  simplified  implementation  of  a  previously 
developed  multiple  source  model  which  can  compress  the 
Monte  Carlo  phase  space  data  by  a  factor  of  1000  or 
more.^^“^^  When  commissioning  another  accelerator  of  the 
same  type,  the  energy  spectra  in  the  beam  model  are  tuned  to 
match  standard  measured  data  such  as  depth  doses  and  dose 
profiles.  Using  this  approach,  we  do  not  have  to  simulate 
every  accelerator  individually.  Only  one  reference  accelera¬ 
tor  needs  to  be  simulated  for  a  type  of  accelerator,  and  this 
can  be  done  carefully  by  someone  with  Monte  Carlo  exper¬ 


tise.  In  this  paper,  a  Varian  Clinac  2100C  accelerator  is  cho¬ 
sen  as  the  reference  machine.  The  machine  is  simulated  us¬ 
ing  the  BEAM  code"^"^  and  a  four-source  beam  model  is 
established  based  on  the  simulated  beam  phase  space  infor¬ 
mation,  The  accuracy  of  the  Monte  Carlo  dose  distributions 
calculated  with  the  model  is  verified.  Then,  the  model  based 
on  the  reference  beam  is  used  to  commission  three  other 
electron  beams.  Two  beams  are  also  from  the  reference  ma¬ 
chine  but  with  incident  energies  significantly  different  from 
that  of  the  reference  beam.  The  third  beam  is  from  another 
Clinac  2100C  machine  at  a  different  institution.^^  The  valid¬ 
ity  of  the  proposed  commissioning  approach  is  demonstrated 
by  commissioning  these  three  beams. 

II.  METHODS  AND  MATERIAL 
A.  Beam  modeling 
1.  General  considerations 

Beam  modeling  is  the  first  step  in  our  hybrid  commission¬ 
ing  procedure  for  a  Monte  Carlo  treatment  planning  system. 
A  beam  model  for  a  type  of  accelerator  is  established  using 
the  Monte  Carlo  simulated  phase  space  information  for  the 
reference  beam.  The  beam  data  are  modeled  using  the  mul¬ 
tiple  source  model  developed  by  Ma  et  which  is 

modified  in  the  current  work  for  use  in  the  commissioning 
procedure.  The  major  modifications  of  the  model  are  dis¬ 
cussed  here. 

The  multiple  source  model  is  based  on  the  observation 
that  particles  from  different  components  of  an  accelerator 
have  significantly  different  energy,  angular,  and  spatial  dis¬ 
tributions,  while  the  particles  from  the  same  component  have 
very  similar  characteristics.^^"^^  Therefore,  the  particles 
from  different  components  of  an  accelerator  can  be  treated  as 
they  are  from  different  sub-sources.  Each  sub-source  repre¬ 
sents  a  critical  component  in  the  treatment  head  and  its  geo¬ 
metrical  dimensions  are  determined  by  the  component  di¬ 
mensions.  Each  sub-source  has  its  own  energy  spectrum  and 
planar  fluence  distribution  derived  from  the  simulated  phase 
space  data.  When  the  model  is  used  for  dose  calculation,  the 
incident  energy  and  position  of  a  particle  are  sampled  from 
the  corresponding  stored  energy  spectrum  and  planar  fluence 
distribution.  The  incident  direction  of  the  particle  is  recon¬ 
structed  by  sampling  the  position  of  the  particle  on  the  sub¬ 
source  and  on  the  phantom  surface.  The  correlation  between 
the  particle  position  and  incident  angle  is  naturally  retained. 

A  primary  reason  to  develop  the  multiple  source  model 
was  to  find  a  concise  way  to  replace  the  huge  phase  space 
data  files  generated  from  Monte  Carlo  simulations.^^’^^  The 
emphasis  of  the  current  work  is  to  develop  a  clinically  prac¬ 
tical  commissioning  procedure  for  Monte  Carlo  treatment 
planning.  The  multiple  source  model  is  simplified  to  make 
the  commissioning  procedure  as  simple  as  possible  while 
trying  to  maintaining  dose  calculation  accuracy  under  all  cir¬ 
cumstances  of  clinical  relevance.  The  number  of  sub-sources 
in  the  model  is  minimized  and  only  those  sub-sources  of 
dosimetric  significance  are  retained.  The  dependence  of  the 
model  on  the  detailed  information  of  accelerators  is  reduced. 
Sub-sources  are  represented  by  dimensionless  geometric  ob- 
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jects,  such  as  points  and  lines,  instead  of  the  actual  geometri¬ 
cal  shapes  and  sizes  of  the  treatment  head  components  as 
used  previously.^^’^^ 

Ma  et  ai  established  their  multiple  source  model  from  the 
Monte  Carlo  simulated  phase  space  data  on  the  patient  sur¬ 
face  (at  100  cm  SSD)  and  the  last  scraper  of  the  electron 
applicator  was  included  in  the  model  as  a  sub-source.^^’^^  In 
this  work,  the  treatment  head  is  simulated  using  the  BEAM 
code  down  to  just  above  the  last  scraper  of  the  electron 
applicator, where  the  patient  specific  cutout  is  inserted.  The 
last  scraper,  as  well  as  the  field-defining  cutout,  are  simu¬ 
lated  together  with  the  patient  CT  phantom  when  performing 
Monte  Carlo  treatment  planning  dose  calculations.  The  ad¬ 
vantage  of  this  method  is  that  the  beam  model  is  patient 
independent.  However,  this  approach,  compared  to  that  of 
Ma  et  requires  more  careful  beam  modeling  since  an 

air  gap  between  the  last  scraper  and  patient  surface  is  not 
included  in  the  original  BEAM  simulation. 

In  the  present  paper,  the  beam  modeling  approach  is  ap¬ 
plied  to  a  Varian  Clinac  2100C  machine  at  Stanford  Medical 
Center,  which  is  chosen  as  the  reference  machine  to  build 
beam  models.  At  first,  a  very  detailed  model  for  each  beam  is 
used  as  the  starting  point.  All  the  critical  components  of  the 
treatment  head  are  modeled  as  sub-sources.  With  this  model, 
the  phase  space  information  of  the  electron  beam  can  be 
precisely  reconstructed  and  the  dose  distribution  in  a  water 
phantom  can  be  accurately  calculated.  Then,  the  number  of 
sub-sources  is  gradually  reduced  while  maintaining  the  accu¬ 
racy  in  dose  distribution  calculation.  We  find  that,  a  point 
electron  source  with  the  energy  spectrum  obtained  from  the 
Monte  Carlo  simulation  is  able  to  give  a  reasonably  accurate 
depth-dose  curve,  which  is  consistent  with  the  previous 
observation.^^  By  adding  another  point  photon  source,  the 
bremsstrahlung  tail  in  the  depth  dose  distribution  can  be  re¬ 
produced  accurately.  However,  it  is  found  that  the  penumbra 
at  the  phantom  surface  generated  with  this  two  point  source 
model  is  sharper  compared  to  that  generated  with  the  full 
phase  space  data.  In  order  to  get  the  dose  profiles  correct,  we 
find  that,  in  addition  to  the  two  point  sources,  two  square 
ring  electron  sources  (which  emit  electrons  isotropically)  are 
needed  to  represent  electrons  scattered  from  the  applicator 
scrapers.  (The  term  square  ring  is  used  here  to  represent  the 
edge  of  a  square.)  Therefore,  the  beam  model  should  include 
four  sub-sources:  a  point  electron  source  for  direct  electrons 
(which  do  not  interact  with  the  beam  defining  system)  and 
electrons  scattered  from  the  primary  collimator,  movable 
jaws  and  shieldings,  a  point  photon  source  for  all  contami¬ 
nant  photons,  and  two  square  ring  electron  sources  for  elec¬ 
trons  scattered  from  the  first  two  scrapers  of  the  electron 
applicator  (the  third  also  last  scraper  is  not  included  in  the 
model). 

As  described  previously,  beam  modeling  consists  of  two 
steps,  namely,  beam  representation  and  beam 
reconstruction. In  beam  representation,  parameters  in  the 
model  are  extracted  from  the  simulated  phase  space  file.  In 
the  current  simplified  model,  these  parameters  include  the 
positions  and  relative  intensities  of  the  sub-sources,  the  en¬ 
ergy  spectra  for  particles  inside  and  outside  the  field  for  each 


sub-source,  and  the  planar  fluence  distributions  on  the  scor¬ 
ing  plane  (directly  above  the  last  scraper)  for  each  sub¬ 
source.  Beam  reconstruction  is  performed  when  using  the 
model  for  dose  calculations.  The  phase  space  information  for 
each  particle,  including  the  energy,  position  and  direction,  is 
reconstructed  from  the  scored  source  parameters. 

2.  Beam  representation 

The  positions  of  the  virtual  electron  and  photon  point 
sources  can  be  determined  using  a  method  described  in  Ref. 
36,  which  is  analogous  to  the  pinhole  method.^^  A  thin  an¬ 
nular  aperture  is  selected  on  the  scoring  plane  and  phase 
space  particles  are  allowed  to  pass  through  the  aperture  and 
form  an  image  at  a  distance  below  the  scoring  plane.  Ray 
lines  drawn  through  the  center  of  the  aperture  and  the  peak 
of  the  aperture  image  form  a  virtual  focal  spot,  which  is 
adopted  as  the  position  of  the  point  source.  This  pinhole 
method  is  very  effective  for  the  photon  point  source.  How¬ 
ever,  we  find  that  for  electrons,  this  method  is  only  appli¬ 
cable  to  high  energy  beams.  For  lower  energy  beams,  e.g.,  6 
MeV,  the  virtual  SSD  determined  with  this  method  is  greatly 
overestimated  and  dependent  on  the  radius  of  the  thin  annu¬ 
lar  aperture  on  the  scoring  plane.  This  is  due  to  the  fact  that 
the  in-air  multiple  scattering  of  low  energy  electrons  is  sig¬ 
nificant.  To  overcome  this  problem  we  performed  another 
Monte  Carlo  simulation  of  the  accelerator  treatment  head  by 
replacing  the  intervening  air  with  vacuum.  Then,  based  on 
the  new  phase  space  data,  this  pinhole  method  can  be  used  to 
generate  the  correct  position  for  the  virtual  electron  point 
source,  which  is  independent  of  the  sampling  radius.  The 
effect  of  in-air  multiple  scattering  is  taken  into  account  dur¬ 
ing  beam  reconstruction  by  adding  a  perturbation  to  the  elec¬ 
tron  incident  direction,  as  discussed  later. 

Two  square  ring  sources  of  electrons  are  located  at  the 
corresponding  positions  of  the  two  applicator  scrapers.  The 
sides  of  the  square  rings  correspond  to  the  actual  openings  of 
the  scrapers. 

The  energy  spectrum  for  each  sub-source  is  derived  from 
the  simulated  phase  space  data.  It  was  found  that  the  mean 
energy  of  the  electrons  is  relatively  uniform  inside  the  field 
as  well  as  outside  the  field.  The  change  of  mean  energy  with 
the  distance  from  the  central  axis  is  more  like  a  step 
f unction. Therefore,  in  the  current  model,  each  sub-source 
has  two  different  energy  spectra,  one  for  electrons  inside  and 
the  other  for  electrons  outside  the  treatment  field.  Parameters 
stored  in  the  model  are  the  minimum  and  maximum  ener¬ 
gies,  number  of  energy  bins  as  well  as  the  relative  fluence  for 
each  energy  bin.  The  minimum  and  maximum  energies  cor¬ 
respond  to  the  cutoff  energies  (ECUT  or  PCUT)  and  the 
incident  energy  used  in  the  accelerator  simulation.  The  num¬ 
ber  of  bins  is  determined  by  the  desired  resolution.  For  ex¬ 
ample,  if  we  want  the  uncertainty  in  the  calculated  depth  of 
50%  dose,  /?5o,  to  be  less  than  1%,  the  uncertainty  in  the 
peak  position  of  the  energy  spectrum  should  be  within  1% 
and  therefore  the  bin  width  should  be  smaller  than  1%  of  the 
peak  energy.  For  the  12  MeV  beam,  we  used  128  bins  and 
then  the  bin  width  is  less  than  0.1  MeV.  This  bin  width  is 
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also  found  to  be  small  enough  to  represent  the  peak  width  of 
the  energy  spectrum,  which  has  a  significant  effect  on  the 
slope  of  the  depth  dose  fall-off. 

The  planar  fluence  distribution  for  each  sub-source  is  also 
derived  from  the  simulated  phase  space  data  and  recorded  on 
the  scoring  plane  using  a  grid  scheme.  Within  each  pixel  of 
the  grid,  the  planar  fluence  is  assumed  to  be  uniform.  The 
dimension  of  a  pixel  is  dependent  on  field  size,  usually  from 
1  to  3  mm.  Parameters  used  to  represent  the  planar  fluence 
distribution  include  the  treatment  field  dimension,  the  num¬ 
ber  of  pixels  and  the  relative  intensity  of  each  pixel,  for  each 
sub-source.  It  has  been  found  that  in  general,  the  mean  en¬ 
ergy  varied  from  position  to  position  in  the  treatment  field  by 
less  than  10%  for  a  given  sub-source.^^  Thus,  it  is  reasonable 
to  store  and  then  sample  the  particle  energy  and  position 
independently. 

The  angular  distributions  are  not  scored  explicitly.  They 
are  reconstructed  during  the  dose  calculation,  as  described  in 
the  next  section. 

Finally,  the  simulated  phase  space  information  is  repre¬ 
sented  with  a  set  of  parameters  for  each  sub- source.  The 
resultant  source  parameter  file  is  much  smaller  (about 
100  kilobytes)  than  the  original  phase  space  file 
(>100  megabytes). 

3.  Beam  reconstruction 

When  performing  dose  calculation  in  a  patient’ s  CT  phan¬ 
tom,  the  source  parameter  file  is  used  to  reconstruct  the 
phase  space  information  (energy,  position  and  direction)  of 
every  incident  particle.  The  beam  reconstruction  process 
consists  of  the  following  steps: 

(1)  Determine  from  which  particular  sub-source  a  particle 
originates  by  sampling  from  the  relative  source  intensity 
of  each  sub-source. 

(2)  Determine  the  position  on  the  sub-source  (excluding 
point  sources)  where  the  particle  is  emitted. 

(3)  Sample  the  particle  position  on  the  scoring  plane  from 
the  fluence  pattern  for  the  sub-source. 

(4)  Sample  the  particle  energy  from  the  energy  spectrum  for 
the  given  sub-source  based  on  the  particle  position  (in¬ 
side  or  outside  the  treatment  field). 

(5)  Determine  the  particle  incident  angle  by  connecting  the 
position  on  the  sub-source  from  where  the  particle  is 
emitted  and  the  position  of  the  particle  on  the  scoring 
plane. 

(6)  Add  in-air  perturbation  on  the  particle  direction  if  it  is  an 
electron. 

The  sampling  from  the  relative  source  intensity  distribu¬ 
tion  is  done  by  a  table-look-up  method.^^  All  the  sub-sources 
are  in  turn  numbered  from  1  to  N  (here  = 4  for  accelerators 
with  designs  similar  to  the  Varian  Clinac  2100C  machine) 
and  the  relative  intensity  of  the  /th  sub-source  is  /?,(/ 
=  1,...,A).  The  accumulative  source  intensity  for  the  /th 
sub-source,  P,(  =  2j^j/?y),  is  multiplied  by  a  large  integer 
M.  The  value  of  M  is  determined  according  to  the  desired 
sampling  precision  of  the  relative  source  intensity.  For  ex- 

Medical  Physics,  Vol.  27,  No.  1,  January  2000 


(-a/2.ay2)  egde2  (a/2, a/2) 


Fig.  1 .  A  diagram  for  illustrating  the  sampling  algorithm  from  a  square  ring 
source.  The  origin  of  the  coordinate  system  is  at  the  center  of  the  square 
ring. 


ample,  M=  1000  corresponds  to  the  precision  of  0.1%  in 
sampling  from  the  relative  source  intensity  distribution. 
Then,  a  one  dimensional  array  of  M  elements  is  prepared  by 
assigning  value  i  to  array  elements  from  /A^r(MP,_,)  to 
INT(MPj),  where  the  operator  I  NT  returns  the  integer  part 
of  a  real  number.  During  the  beam  reconstruction,  a  random 
integer  number  K  between  1  and  M  is  generated  and  the 
value  of  the  A^th  array  element  is  the  sub-source  number 
where  the  particle  is  emitted.  Such  a  table-look-up  method  is 
of  very  high  sampling  speed  and  efficiency.  Its  speed  is  also 
independent  of  the  number  of  sub-sources.  The  sampling 
precision  is  usually  adequate  as  long  as  an  large  enough  ar¬ 
ray  is  used. 

According  to  the  energy  spectra  on  the  scoring  plane,  the 
bin  number  in  which  the  particle  energy  falls  is  also  sampled 
using  the  table-look-up  method.  An  additional  uniform  sam¬ 
pling  is  done  within  the  given  energy  bin  to  make  the  particle 
energy  continuous. 

The  same  table-look-up  method  is  also  used  to  sample  the 
pixel  number  corresponding  to  a  particle  position  on  the 
scoring  plane.  Another  sampling  is  performed  uniformly  to 
determine  the  particle’s  coordinates  within  the  chosen  pixel. 

For  the  point  sources,  the  incident  angle  of  the  particle  is 
determined  by  constructing  a  ray  line  from  the  point  source 
to  the  position  of  the  particle  on  the  scoring  plane.  For  the 
square  ring  electron  sources,  we  need  to  determine  where  the 
electron  comes  from  on  the  ring.  This  is  done  in  two  steps. 
First,  we  determine  from  which  edge  of  the  square  ring  the 
electron  is  emitted.  Second,  we  determine  from  where  on  the 
chosen  edge  the  electron  comes. 

As  illustrated  in  Fig.  1,  a  square  ring  of  size  aXa  is 
located  above  the  scoring  plane  at  a  distance  d.  We  assume 
that  the  square  ring  edge  emits  electrons  uniformly  and  iso¬ 
tropically.  Under  this  approximation,  the  probability  for  an 
electron  on  the  scoring  plane  to  come  from  a  point  on  the 
ring  is  proportional  to  the  inverse  square  of  the  distance  be- 
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tween  the  points  on  the  ring  and  plane.  This  approximation 
greatly  simplifies  the  sampling  process.  It  is  found  that  the 
angular  distribution  of  the  scattered  electrons  reconstructed 
based  on  the  fluence  distribution  on  the  scoring  plane  and  the 
emitting  position  on  the  square  ring  source  is  reasonably 
accurate,  although  electrons  scattered  from  the  applicator 
scraper  mainly  correspond  to  electrons  incident  on  the  verti- 
cal  face  of  the  scraper  and  are  dominantly  forward  directed 
Furthermore,  the  dose  profile  at  the  patient  surface  is  greatly 
influenced  by  electrons  scattered  from  the  last  scraper  or  cut¬ 
outs,  which  are  not  included  in  the  source  model  but  will  be 
accurately  simulated  with  the  patient  CT  phantom. 

For  an  electron  at  position  (jco,yo)  the  scoring  plane 
(see  Fig.  1),  the  probability  for  it  to  have  come  from  the  iih 
edge  is  given  as 
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C4=yl{xo-a/2f  +  d^. 

(6) 

Using  the  probabilities  given  above,  the  edge  from  which  the 
electron  has  come  can  be  sampled.  Then,  the  position  (x,y) 
on  the  chosen  edge  is  further  sampled.  For  edges  1  and  2, 
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where  ^  is  a  random  number  uniformly  distributed  from  0 
to  1. 

After  the  electron  position  on  the  square  ring  is  deter¬ 
mined,  the  connection  of  this  position  to  the  position  on  the 
scoring  plane  gives  the  electron’s  incident  direction,  which 
needs  to  be  additionally  perturbed  to  address  the  in-air  mul¬ 
tiple  scattering. 

In  a  previous  implementation  of  this  model,  the  effect  of 
electron  multiple  scattering  in  air  as  well  as  other  materials 
on  its  path  to  the  scoring  plane  was  taken  into  account  by 
sampling  the  electron  perturbation  angle  from  a  Monte  Carlo 
simulated  angular  distribution.^^  This  angular  distribution 
was  stored  while  performing  the  Monte  Carlo  simulation  for 
the  accelerator  and  only  included  electrons  falling  into  a 
small  region  (e.g.,  of  1  cm  radius)  around  beam  central  axis. 
In  the  current  work,  the  effect  of  in-air  multiple  scattering  is 
considered  more  accurately  using  the  Fermi -Eyges  theory 
The  effect  of  other  materials  is  considered  by  adding  a  pa¬ 
rameter  to  the  standard  deviation  of  the  angular  distribution. 
The  Fermi-Eyges  theory  is  a  well-known  small-angle  theory 
and  can  predict  the  multiple  scattering  effect  of  megavoltage 
electrons  in  air  or  other  heavier  materials  as  long  as  the  elec¬ 
tron  effective  pathlengths  are  small.'^^’'^^ 

Assume  that  an  electron  initially  travels  along  the  z  axis. 
According  to  the  Fermi-Eyges  theory,  the  distributions  for 
the  projections  of  the  polar  angle,  ft  on  the  x-z  plane,  6^ 
and  on  the  y  -  z  plane,  Oy ,  are  both  Gaussian  after  electrons 
travel  a  distance,  and  are  given  as^’^^ 
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where  and  are  the  standard  deviations  for  each 
Gaussian  distribution,  respectively.  In  a  homogeneous  mate¬ 
rial,  and  without  the  presence  of  an  electromagnetic  field, 
both  standard  deviations  should  be  the  same,  so  we  let  a 
=^0-0  =0-0^,  Under  the  small-angle  approximation 

(13) 


therefore  the  polar  angle  6  obeys  a  radial  Gaussian  distribu¬ 
tion  while  the  azimuthal  angle  <f)  is  uniformly  distributed  in 
[0,2  tt].  Hence  the  sampling  method  for  these  two  angles  is 
given  as  follows: 


0=oryl-2  In^i, 


(14) 


0  =  27r^2.  (15) 

where  and  ^2  random  numbers  uniformly  distributed 
in  [0,1]. 

According  to  the  Fermi-Eyges  theory,  a  can  be  calcu¬ 
lated  as^’"^^ 

a^=Ao-A]/A2,  (16) 

where 


Medical  Physics,  Vol.  27,  No.  1,  January  2000 


185 


Jiang,  Kapur,  and  Ma:  Electron  beam  modeling  and  commissioning 


185 


1 

A--  K{l-tydt,  /  =  0,1,2.  (17) 

^  Jo 

Here,  K  is  the  electron  linear  scattering  power  and  /  is  the 
distance  at  which  electrons  travel.  The  electron  linear  scat¬ 
tering  power  can  be  fitted  well  using  a  simple  formula  pro¬ 
posed  by  Werner  et  al 

K{E)  =  aE-^,  (18) 

Using  this  formula  we  fitted  the  linear  scattering  power  data 
in  air  supplied  by  ICRU  Report  35^  and  found  that  a 
=  3,329 X  10~^  rad^/cm  and  1.638.  E  is  the  electron  en¬ 
ergy  in  MeV  and  sampled  from  the  energy  spectrum  at  the 
scoring  plane.  The  energy  loss  of  electrons  in  air  is  usually 
very  small  and  can  be  ignored  when  they  travel  from  the 
virtual  source  to  the  scoring  plane.  The  mean  energy  loss  of 
6  MeV  electrons  after  traveling  100  cm  in  air  is  about  4%  of 
its  initial  energy  (estimated  using  the  stopping  power)  and  it 
is  about  2%  for  20  MeV  electrons.  Therefore  a  can  be  given 
as 

(19) 

which  is  a  function  of  electron  energy  and  the  distance  be¬ 
tween  the  virtual  source  and  the  position  on  scoring  plane. 
During  beam  reconstruction,  according  to  the  sampled  elec¬ 
tron  energy,  positions  on  the  scoring  plane  and  on  the  virtual 
source,  a  can  be  calculated.  Then  using  Eqs.  (14)  and  (15)  6 
and  (f)  are  sampled  and  a  perturbation  is  added  to  the  elec¬ 
tron’s  incident  direction. 

The  perturbation  caused  by  in-air  multiple  scattering  can 
be  directly  calculated  using  Eq.  (19)  for  electrons  from  the 
squaring  ring  sources.  For  direct  electrons,  there  are  other 
accelerator  components  in  their  paths  from  the  virtual  point 
source  to  the  scoring  plane  in  addition  to  the  intervening  air, 
such  as  the  exit  window,  scattering  foil,  monitor  chamber, 
mirror  and  protection  window.  The  effect  of  these  materials 
on  electron  angular  distribution  has  been  mainly  included  in 
the  determination  of  the  virtual  electron  point  source  posi¬ 
tion.  We  also  need  to  take  into  account  the  angular  perturba¬ 
tion  caused  by  these  materials.  If  we  know  precisely  the 
material  and  thickness  of  these  parts,  we  can  calculate  their 
effect  on  o',  as  done  by  Keall  and  Hoban."^*  However,  it  is 
usually  difficult  for  users  to  know  this  information  about 
their  accelerator  when  commissioning  a  Monte  Carlo  treat¬ 
ment  planning  system.  Therefore,  we  introduce  a  factor  k  to 
account  for  the  effect  of  these  materials.  For  direct  electrons, 
a  is  then  given  as 

a'^=\K{E)lk.  (20) 

The  factor  k  is  determined  by  fitting  the  angular  distribution 
calculated  using  Fermi-Eyes  theory  to  that  simulated  with 
the  Monte  Carlo  method  for  direct  electrons.  The  introduc¬ 
tion  of  k  factor  provides  a  potentially  tunable  parameter  in 
the  source  model. 

4.  Model  verification 

The  four-source  model  is  verified  dosimetrically  by  com¬ 
paring  the  dose  distributions  in  a  water  phantom  calculated 


using  the  model  with  those  calculated  using  the  full  phase 
space  data.  Dose  distributions  are  calculated  for  various 
combinations  of  three  electron  energies  (6,  12,  and  20  MeV), 
three  applicator  sizes  (6X6,  10X10,  and  20X20cm^),  and 
two  SSDs  (100  cm  and  120  cm). 

The  measurement  of  electron  beam  applicator  factors  (de¬ 
fined  as  the  ratio  of  the  open  field  dose  in  water  at  for  a 
given  applicator  to  that  of  the  reference  applicator,  typically 
the  10X10  or  15  Xl5cm^,  for  the  same  beam  energy)  is 
done  during  accelerator  commissioning  for  all  energy/cone 
combinations.  Therefore,  the  applicator  factors  will  be  sup¬ 
plied  by  the  user  when  performing  the  model  commission¬ 
ing.  Cutout  factors  (defined  as  the  ratio  of  the  dose  in  water 
at  for  a  blocked  field  to  that  of  the  open  field  for  the 
same  applicator  and  beam  energy)  are  patient  specific  and 
not  always  easy  to  measure  accurately  for  all  clinical  situa¬ 
tions.  Therefore,  the  model  should  be  able  to  calculate  cutout 
factors.  To  demonstrate  this,  we  compare  the  model  calcu¬ 
lated  cutout  factors  with  those  measured  and  calculated  by 
Kapur  et  al  using  a  full  Monte  Carlo  simulation. 


B.  Beam  commissioning 

The  four-source  model  which  is  built  using  a  Varian  Cli- 
nac  2100C  accelerator  as  the  reference  machine  can  be  used 
to  reconstruct  electron  beams  from  other  Clinac  2100C  ac¬ 
celerators  by  tuning  the  energy  spectra  in  the  model. 

For  accelerators  with  exactly  the  same  design,  the  major 
different  is  the  electron  incident  energy  due  to  the  on-site 
tuning  to  suit  the  user.  This  energy  approximately  corre¬ 
sponds  to  the  maximum  energy  of  all  the  stored  energy  spec¬ 
tra  in  the  source  model.  It  is  found  that  the  energy  spectra  for 
all  sub-sources  are  very  similar  for  different  accelerators  of 
the  same  type.  When  the  incident  energy  is  changed,  the 
energy  spectra  can  be  approximated  as  stretched  or  com¬ 
pressed  along  the  energy  axis  accordingly.  The  depth  dose 
curve  is  very  sensitive  to  the  electron  incident  energy  and 
therefore  used  to  adjust  the  maximum  energy,  ihe 

model.  The  relationship  between  the  incident  energy, 
and  /?5o  has  been  studied  by  simulating  the  reference  accel¬ 
erator  using  a  lOX  10  cm^  cone  and  100  cm  SSD  with  vari¬ 
ous  incident  energies.  Then,  the  variation  of  as  a  function 
of  the  variation  of  Rso  is  established  for  this  type  of  accel¬ 
erator.  This  relationship  is  used  as  a  guide  to  tune  the  maxi¬ 
mum  energy  in  the  model  to  commission  a  clinical  beam. 

The  proposed  commissioning  approach  can  be  summa¬ 
rized  as  follows: 

(1)  Chose  an  accelerator  as  the  reference  machine  for  all 
other  accelerators  of  the  same  design,  and  carefully  per¬ 
form  full  Monte  Carlo  simulations  for  the  electron 
beams  of  various  nominal  energies  from  the  reference 
machine  with  lOX  10 cm^  applicator. 

(2)  Build  the  source  models  for  the  simulated  beams  based 
on  the  Monte  Carlo  phase  space  data,  perform  Monte 
Carlo  dose  calculation  in  water  for  100  cm  SSD,  and 
record  the  maximum  energy,  ,  in  the  model  and  the 
/?50^  value  for  each  beam. 
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Fig.  2.  Effects  of  intervening  air  and  sampling  radius  on  virtual  SSD  (de¬ 
fined  here  as  the  distance  from  the  point  source  to  the  scoring  plane)  for  the 
electron  point  source  determined  with  the  pinhole  method,  (a)  6  MeV;  (b)  20 
MeV. 


(3)  For  the  beam  to  be  commissioned,  find  the  value 
of  the  measured  depth-dose  curve  in  water  for  10 
X  lOcm^  applicator  and  100  cm  SSD. 

(4)  Select  the  reference  beam  which  has  the  same  or  closest 
nominal  energy  as  the  commissioning  beam.  Let  /  =  0, 

r(0  D(0_D(re0 

^inax*~^max>  ™  ^5Q”^50  * 

(5)  Calculate  =  If  where  6  is 

the  pre-set  convergence  tolerance,  stop  iteration  and  use 
^max  ^  maximum  energy  in  the  source  model  for  the 
commissioning  beam;  otherwise,  go  to  the  next  step. 

(6)  According  to  the  relationship  between  and  A/?5o» 

calculate  using  A/?g(]  and  then  calculate 

=£®^-A£<'\ 

(7)  Calculate  the  dose  distribution  using  the  source  model 

with  £'mtx^  the  corresponding  . 

(8)  Let  /<—/+!;  go  back  to  step  5. 


clinical  situations.  The  iteration  process  converges  very  fast; 
usually  only  two  or  three  iterations  are  needed  even  for  € 
much  smaller  than  1  mm. 

The  commissioning  approach  has  been  applied  to  three 
electron  beams.  A,  B,  and  C.  The  reference  beam  is  the  same 
for  all  three  beams,  which  is  from  the  reference  Clinac 
2100C  accelerator  with  12.0  MeV.  Beam  A  and  beam 
B  are  also  from  the  reference  machine  but  with  Ei^  as  9.0 
MeV  and  15.0  MeV,  respectively.  These  two  beams  are  used 
to  mimic  two  clinical  beams  of  the  same  nominal  energy  as 
the  reference  beam  but  with  significantly  different  incident 
energies.  Of  course,  in  reality,  the  electron  incident  energy 
will  not  be  tuned  so  much  (±3  MeV).  These  two  beams  are 
used  as  extreme  cases  to  test  the  commissioning  approach. 
Beam  C  is  a  9  MeV  electron  beam  from  another  Clinac 
2100C  accelerator.  The  dose  distributions  for  beam  C  are 
taken  from  the  published  data.^^ 


III.  RESULTS  AND  DISCUSSION 

Figure  2  shows  the  effects  of  intervening  air  and  sampling 
radius  on  the  electron  and  photon  virtual  SSD  determined 
with  the  pinhole  method  for  6  MeV  and  20  MeV  beams.  It 
can  be  seen  that  for  photons  and  high  energy  electrons  (20 
MeV)  effects  of  intervening  air  and  sampling  radius  on  the 
positions  of  virtual  point  sources  are  negligible.  However, 
for  low  energy  electrons  (6  MeV),  these  effects  are  signifi¬ 
cant.  Therefore,  to  obtain  the  accurate  virtual  electron  point 
source  position  for  low  energy  beams,  the  phase  space  simu¬ 
lated  without  intervening  air  should  be  used. 

Figure  3  shows  the  comparison  between  the  angular  dis¬ 
tributions  for  direct  electrons  in  6,  12  and  20  MeV  beams 
calculated  with  the  Fermi -Eyges  theory  and  the  Monte  Carlo 
method.  We  can  see  that,  in  general,  the  fitted  angular  distri¬ 
butions  based  on  the  Fermi -Eyges  theory  match  well  with 
those  calculated  with  the  Monte  Carlo  simulation.  We  also 
notice  that  at  large  angles  the  Fermi -Eyges  theory  slightly 


The  first  two  steps  only  need  to  be  do  once  for  all  accelera¬ 
tors  of  the  same  design.  The  convergence  tolerance,  6,  is  set 
by  the  user,  usually  according  to  the  estimated  measurement 
error  in  R^q  .  For  example,  €=  1  mm  is  good  enough  in  most 


Fig.  3.  Angular  distributions  for  direct  electrons  calculated  using  the 
Fermi-Eyges  theory  and  the  Monte  Carlo  method.  Beam  energies  are  6 
MeV,  12  MeV  and  20  MeV.  The  fitted  k  factor  is  1.540  for  6  MeV,  1.501 
for  12  MeV  and  1.571  for  20  MeV.  Each  distribution  is  normalized  to  have 
unit  area  under  the  curve. 
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Fig.  4.  Dose  distributions  in  water  for  12  MeV  electron  beam  with  10 
XlOcm^  applicator  at  100  cm  SSD,  calculated  with  full  phase  space  data 
and  source  model:  (a)  Depth-dose  distributions;  (b)  dose  profiles  at  depths  of 
2  cm  and  5  cm.  Curves  are  normalized  to  the  dose  at  • 


underestimates  electron  fluence  due  to  the  fact  that  it  is  a 
small-angle  theory.  We  found  that  the  small  discrepancy 
does  not  have  any  significant  effect  on  the  final  dose  distri¬ 
butions.  Therefore,  the  Fermi-Eyges  theory  with  the  fitted  k 
factor  can  be  used  to  account  for  the  angular  perturbations  of 
electrons  on  their  way  from  the  source  to  the  scoring  plane. 

The  four-source  model  was  tested  by  comparing  the  dose 
distributions  calculated  by  the  model  with  those  calculated 
by  full  phase  space  data  for  various  combinations  of  three 
electron  energies  (6,  12,  and  20  MeV),  three  applicator  sizes 
(6X6,  10X10,  and  20X20cm^),  and  two  SSDs  (100  cm  and 
120  cm).  For  all  the  cases  tested,  the  agreement  of  l%-2%f 
1-2  mm  has  been  achieved.  Figure  4  shows  the  comparison 
for  a  12  MeV  beam  with  a  lOX  10  cm^  cone  at  100  cm  SSD. 
Figure  5  gives  the  depth-dose  curves  and  dose  profiles  for 
20  MeV  beam  with  6X6 cm^  cone  at  120  cm  SSD,  calcu¬ 
lated  with  both  the  source  model  and  full  phase  space  data. 
All  the  curves  in  Figs.  4  and  5  are  normalized  to  the  doses  at 
^max  •  The  Monte  Carlo  uncertainty  is  always  less  than  0.5% 
and  therefore  not  shown  on  the  curves.  In  both  figures  the 
agreement  between  the  full  Monte  Carlo  simulations  and  the 
source  model  calculations  is  better  than  1%/1  mm.  Keep  in 


Fig.  5.  Dose  distributions  in  water  for  20  MeV  , electron  beam  with  6 
X6  cm^  applicator  at  120  cm  SSD,  calculated  with  full  phase  space  data  and 
source  model:  (a)  Depth-dose  distributions;  (b)  dose  profiles  at  depths  of  2 
cm  and  7  cm.  Curves  are  normalized  to  the  dose  at  ♦ 


mind  that  20  cm  air  gap  is  rarely  used  in  actual  clinical 
situations.  Therefore,  we  have  demonstrated  that  the  simpli¬ 
fied  four-source  model  can  be  used  for  accurate  dose  calcu¬ 
lations  even  for  extreme  cases  (such  as  very  large  extended 
SSDs). 

The  capability  of  the  model  for  calculating  the  relative 
beam  output  was  also  tested.  Table  I  shows  cutout  factors  for 
various  square  inserts  in  a  lOXlOcm^  applicator  for  6,  12 
and  20  MeV  electron  beams.  It  is  found  that  the  cutout  fac¬ 
tors  calculated  with  the  four-source  model  are  within  about 
±2%  compared  to  the  measured  values  except  for  one  case 
where  we  see  2.5%  difference.  This  is  at  about  the  same 
accuracy  level  as  the  full  Monte  Carlo  simulation  and  con¬ 
sidered  to  be  acceptable  for  clinical  use. 

The  relationship  between  and  R^q  for  the  reference 
accelerator  with  a  lOXlOcm^  cone  and  100  cm  SSD  is 
shown  in  Fig,  6,  A  linear  relationship  was  found  and  fitted  as 

i?in=2.597/?5o+ 0,633,  (21) 

It  gives  the  relationship  between  the  variation  in  the  incident 
energy  and  the  variation  in  ^50  as 

2.597 A/?5o.  (22) 
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Table  I.  The  electron  cutout  factors  for  various  square  inserts  in  10X10  cm^  applicator  for  6,  12,  and  20  MeV 
beams  calculated  with  the  source  model  and  the  full  phase  space  data  and  compared  with  the  measurement  of 
Kapur  et  al  (Ref.  45).  The  values  in  parenthesis  indicate  the  difference  of  the  data  calculated  with  the  source 
model  or  the  full  phase  space  from  the  measured  data. 


Energy 

(MeV) 

Insert 

(cm^) 

Cutout  factor 

Source  model 

Full  phase  space 

Measurement 

6 

2X2 

0.803  (2.5%) 

0.765  (-1.3%) 

0.778 

3X3 

0.930  (0.3%) 

0.923  (-0.4%) 

0.927 

4X4 

0.970  (-1.8%) 

0.982  (-0.6%) 

0.988 

8X8 

1.002  (-0.1%) 

1.005  (-0.2%) 

1.003 

12 

2X2 

0.881  (-0.8%) 

0.861  (-2.8%) 

0.889 

3X3 

0.908  (-2.0%) 

0.930  (0.2%) 

0.928 

4X4 

0.942  (-2.1%) 

0.956  (-0.7%) 

0.963 

8X8 

0.999  (0.8%) 

1.002  (1.1%) 

0.991 

20 

2X2 

0.963  (-1.1%) 

0.957  (-1.9%) 

0.976 

3X3 

0.989  (-0.4%) 

0.968  (-2.5%) 

0.993 

4X4 

0.993  (-1.8%) 

0.993  (-1.8%) 

1.011 

8X8 

0.999  (-0.5%) 

0.993  (-1.1%) 

1.004 

Equation  (22)  is  used  for  tuning  the  maximum  energy  in  the 
source  model  to  match  the  measured  depth  dose  curves  when 
commissioning  a  clinical  beam. 

Figure  7  shows  the  dose  distributions  for  the  reference 
beam,  beam  A,  and  beam  B  with  the  applicator  size  of  10 
X  10  cm^  and  SSD  of  100  cm.  All  the  curves  are  normalized 
to  the  dose  at  .  The  statistical  uncertainty  (Icr)  in  all  the 
Monte  Carlo  dose  calculations  was  kept  to  be  smaller  than 
0.5%,  therefore,  the  error  bars  are  smaller  than  the  symbol 
size  and  not  shown  on  the  curves.  The  maximum  energy  in 
the  source  model  was  adjusted  to  8.87  MeV  to  match  the 
dose  distributions  of  the  beam  A  (£;,„=  9.0  MeV)  and  to 
15,17  MeV  to  match  the  dose  distributions  of  the  beam  B 
15.0 MeV).  The  difference  between  the  depth-dose 
curves  calculated  by  the  adjusted  models  and  the  full  Monte 
Carlo  simulation  is  always  less  than  0.5%  for  both  beam  A 
and  beam  B.  For  dose  profiles,  the  difference  is  usually  less 
than  1%  except  that  in  the  shoulder  region  for  beam  B  the 
difference  is  about  2%. 

Figure  8  shows  the  dose  distributions  for  the  reference 
beam  and  beam  C  with  the  applicator  size  of  1  OX  10  cm^  and 
SSD  of  100  cm.  Again,  the  curves  are  normalized  to  the  dose 
at  ^/max  arid  the  Monte  Carlo  uncertainty  is  lower  than  0,5%. 
In  this  case,  the  maximum  energy  in  the  source  model  was 
adjusted  to  1 1 .25  MeV.  The  dose  distributions  calculated  by 
the  source  model  with  the  adjusted  maximum  energy  agree 
very  well  (1%/1  mm)  with  the  published  data.^^ 

Table  II  gives  and  for  the  reference  and  Monte 
Carlo  simulated  beams,  and  and  /?5o  for  the  adjusted 
source  models.  For  the  reference  beam,  E^  was  directly 
obtained  from  the  full  Monte  Carlo  simulation.  For  beam  C, 
E'in  is  unknown.  In  this  study,  we  set  6=0.01  cm.  Therefore, 
the  /?5o’s  calculated  using  the  adjusted  source  model  match 
with  the  full  Monte  Carlo  simulation  to  within  0.01  cm.  Of 
course,  we  will  not  use  such  a  small  e  in  real  clinical  appli¬ 
cations  since  it  is  much  smaller  than  the  measurement  uncer¬ 
tainty  in  /?5o.  Here,  we  just  want  to  demonstrate  the  capa¬ 


bility  of  the  method  to  reproduce  accurately. 

We  have  applied  the  commissioning  approach  to  electron 
beams  from  a  Clinac  2300C/D  accelerator  in  our  institution. 
The  reference  machine  is  still  the  same  Clinac  2100C  accel¬ 
erator.  These  two  machines  are  sufficiently  similar  to  each 
other  in  treatment  head  geometry.  Their  dosimetric  charac¬ 
teristics  are  very  close  to  each  other  due  to  the  beam  tuning 
during  linac  acceptance.  Therefore,  it  is  not  surprising  to  see 
that  the  dose  distributions  calculated  with  the  adjusted  source 
model  agree  well  (within  l%-2%  or  1-2  mm)  with  the  mea¬ 
sured  data. 

These  preliminary  results  have  shown  that  the  proposed 
hybrid  commissioning  approach  can  be  used  for  accelerators 
of  the  same  design  to  account  for  the  dosimetric  variations 
mainly  caused  by  the  on-site  tuning  of  the  electron  incident 
energy.  The  capability  of  the  approach  to  handle  large  varia¬ 
tion  in  the  electron  incident  energy  has  been  demonstrated.  It 
is  believed  that  for  most  clinical  accelerators  of  the  same 
type,  their  treatment  head  designs  are  exactly  the  same  or  at 
least  very  similar,  therefore  the  dosimetric  difference  can 
usually  be  traced  back  to  the  difference  in  the  electron  inci¬ 
dent  energy.  Therefore,  the  current  approach  should  be  ap¬ 
plicable  in  most  clinical  situations.  In  the  future  work,  the 
method  will  be  evaluated  under  more  critical  conditions, 
such  as  small  field  sizes,  extended  SSD,  and  heterogeneous 
phantoms. 

The  general  idea  proposed  here  should  also  work  for  other 
types  of  accelerators,  although  we  have  selected  the  Varian 
Clinac  2100C  accelerators  in  the  current  study.  For  each  type 
of  accelerator,  a  reference  machine  should  be  carefully  simu¬ 
lated  using  the  Monte  Carlo  method.  A  source  model,  which 
may  consist  of  a  different  number  of  sub-sources,  can  be 
established  based  on  the  simulated  data.  Then,  the  maximum 
energy  in  the  model  can  be  adjusted  to  commission  electron 
beams  from  other  accelerators  of  the  same  type. 

In  some  situations,  the  proposed  commissioning  approach 
may  not  be  directly  applicable.  For  example,  the  measured 
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Electron  mean  range  (cm) 


Fig.  6.  The  relationship  between  electron  incident  energy  at  accelerator  exit 
window  (Ejn)  and  the  depth  of  50%  dose  (R^o)  for  electron  beams  from  the 
reference  accelerator  with  10X10  cm^  applicator  and  100  cm  SSD.  Symbols 
are  calculated  with  Monte  Carlo  simulation  of  the  accelerator  treatment 
head.  Solid  line  is  the  fitted  result  with  formula  £'in=  2.597/? 50+ 0.633. 


dose  distributions  used  for  commissioning  more  or  less  con¬ 
tain  measurement  errors,  depending  on  the  measurement 
techniques  and  the  experience  of  the  person  who  performs 
the  measurements.  Since  only  the  maximum  energy  is  the 
adjustable  parameter  in  the  current  source  model,  our  ap¬ 
proach  may  not  be  able  to  exactly  match  the  measured  data. 
Occasionally,  an  accelerator  used  in  the  clinic  may  differ 
from  its  original  design  in  addition  to  the  electron  incident 
energy.  Some  parts  in  the  accelerator  treatment  head  may  be 
replaced  with  nonstandard  ones.  In  this  case,  we  can  always 
perform  a  full  Monte  Carlo  simulation  for  this  unique  accel¬ 
erator  and  build  its  own  source  model.  Alternatively,  we  can 
make  the  present  approach  more  versatile  to  handle  those 
situations.  More  parameters  in  the  source  model  other  than 
the  maximum  energy,  such  as  the  relative  intensity  of  each 
sub-source,  the  k  factor  for  in-air  perturbation  for  the  direct 
electrons,  and  the  field  size,  can  be  adjusted  to  match  the 
measured  dose  distributions.  For  example,  the  adjustment  of 
the  relative  intensity  of  the  photon  source  will  ensure  a  good 
match  to  the  bremsstrahlung  tail  in  the  depth-dose  curve.  If 
some  of  the  materials  in  the  paths  of  direct  electrons,  such  as 
the  scattering  foil,  monitor  chamber  or  mirror,  are  different 
from  those  used  in  the  reference  accelerator,  the  adjustment 
of  the  k  factor  can  yield  a  better  estimation  of  the  electron 
angular  perturbation.  The  adjustment  of  the  field  size  in  the 
source  model  should  recover  the  measurement  error  in  the 
width  of  the  dose  profiles  (e.g.,  errors  of  the  order  of  about  1 
mm  are  not  rare  in  a  clinical  situation).  In  summary,  the 
introduction  of  more  adjustable  parameters  in  the  source 
model  will  make  the  current  commissioning  approach  more 
powerful.  This  possibility  will  be  investigated  in  our  future 
study. 

IV.  CONCLUSIONS 

A  hybrid  commissioning  approach  based  on  a  multiple 
source  model  has  been  proposed  for  Monte  Carlo  treatment 
planning.  It  has  been  demonstrated  that  a  simplified  four- 
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Fig.  7.  Dose  distributions  in  water  for  electron  beams  from  the  reference 
accelerator  with  lOX  10  cm^  applicator  and  at  100  cm  SSD.  The  reference 
beam  has  the  electron  incident  energy  of  12.0  MeV.  A  source  model  was 
built  based  on  the  Monte  Carlo  simulation  of  the  reference  beam.  The  maxi¬ 
mum  energy  of  the  energy  spectra  in  the  model  was  adjusted  to  8.87  MeV  to 
match  the  beam  A  (with  £:in=9.0MeV)  and  15.17  MeV  to  match  the  beam 
B  (with  15.0  MeV).  Lines  are  dose  distributions  from  the  full  Monte 
Carlo  simulations.  Symbols  are  data  calculated  by  the  source  model  with 
adjusted  maximum  energies.  All  data  are  normalized  to  the  doses  at  • 
(a)  Depth-dose  distributions;  (b)  dose  profiles  at  depths  of  2  cm  and  3  cm  for 
beam  A;  (c)  dose  profiles  at  depths  of  3  cm  and  5  cm  for  beam  B. 


source  model  can  be  used  to  generate  accurate  Monte  Carlo 
dose  distributions  for  electron  beams  from  Varian  Clinac 
2100c  accelerators.  The  model  includes  a  point  electron 
source  for  direct  electrons  and  electrons  scattered  from  pri- 
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Fig.  8.  Dose  distributions  in  water  for  electron  beams  with  lOxiOcm^ 
applicator  and  at  100  cm  SSD.  The  reference  beam  (solid  line)  is  from  a 
Clinac  2100C  accelerator  with  Ein-12.0MeV.  The  dose  distributions  for 
beam  C  (dashed  lines)  is  taken  from  the  published  data  for  a  9  MeV  beam 
from  another  Clinac  2100C  accelerator  with  type  III  applicator  (Ref.  38).  A 
source  model  built  based  on  the  Monte  Carlo  simulation  of  the  reference 
beam  is  used  to  match  beam  C  by  adjusting  the  maximum  energy  to  1 1 .25 
MeV  (open  circles).  All  data  are  normalized  to  the  dose  at  .  (a)  Depth- 
dose  distributions;  (b)  dose  profiles  at  depth  of  0.1  cm. 

mary  collimator  and  jaws,  a  point  photon  source  for  all  con¬ 
taminant  bremsstrahlung  photons,  and  two  square  ring  elec¬ 
tron  sources  representing  electrons  scattered  from  two  scrap¬ 
ers  (other  than  the  last  scraper)  of  the  Varian  electron 
applicator  (type  III).  It  was  found  that  the  position  of  the 
virtual  point  source  can  be  determined  accurately  using  the 
pinhole  method  for  photons  in  all  cases  and  electrons  in  high 
energy  beams.  For  low  energy  beams,  we  should  use  the 
Monte  Carlo  phase  space  data  which  are  obtained  by  simu¬ 
lating  the  treatment  head  with  the  intervening  air  replaced  by 
vacuum.  We  also  found  that  the  in-air  perturbation  on  the 
electron  incident  direction  can  be  properly  accounted  for  us¬ 
ing  the  Fermi-Eyges  model.  The  source  model  which  was 
built  based  on  the  simulated  phase  space  data  for  the  refer¬ 
ence  accelerator  can  be  used  for  other  accelerators  of  the 
same  type,  by  simply  adjusting  the  energy  spectra  in  the 
model.  The  capability  of  this  commissioning  approach  for 
handling  large  variation  in  the  electron  incident  energy  has 
been  demonstrated. 


Table  II.  Some  parameters  for  the  full  Monte  Carlo  simulations  and  the 
adjusted  source  models.  is  the  electron  incident  energy  at  the  exit  win¬ 
dow  in  the  simulation.  is  the  maximum  energy  in  the  source  model.  R50 
is  the  depth  of  50%  dose  of  the  depth  dose  curve  in  water  for  lOX  10  cm^ 
applicator  and  100  cm  SSD.  A  source  model  was  built  based  on  the  full 
Monte  Carlo  simulation  of  the  reference  beam.  For  beams  to  be  commis¬ 
sioned  (A,  B,  and  C),  £niax  tit  the  model  was  adjusted  to  match  the  corre¬ 
sponding  /?5o  from  the  Monte  Carlo  simulation. 


Beam 

MC  simulation 

Source  model 

tested 

(MeV) 

/?50  (cm) 

£ma,  (MeV) 

/?50  (cm) 

Reference 

12.0 

4.397 

12.00 

4.397 

Beam  A 

9.0 

3.197 

8.87 

3.192 

Beam  B 

15.0 

5.622 

15.17 

5.624 

Beam  C 

unknown 

4.116 

11.25 

4.119 
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Abstract.  The  purpose  of  this  work  was  to  use  Monte  Carlo  simulations  to  verify  the  accuracy  of 
the  dose  distributions  from  a  commercial  treatment  planning  optimization  system  (Corvus,  Nomos 
Corp.,  Sewickley,  PA)  for  intensity-modulated  radiotherapy  (IMRT).  A  Monte  Carlo  treatment 
planning  system  has  been  implemented  clinically  to  improve  and  verify  the  accuracy  of  radiotherapy 
dose  calculations.  Further  modifications  to  the  system  were  made  to  compute  the  dose  in  a  patient  for 
multiple  fixed-gantry  IMRT  fields.  The  dose  distributions  in  the  experimental  phantoms  and  in  the 
patients  were  calculated  and  used  to  verify  the  optimized  treatment  plans  generated  by  the  Corvus 
system.  The  Monte  Carlo  calculated  IMRT  dose  distributions  agreed  with  the  measurements  to 
within  2%  of  the  maximum  dose  for  all  the  beam  energies  and  field  sizes  for  both  the  homogeneous 
and  heterogeneous  phantoms.  The  dose  distributions  predicted  by  the  Corvus  system,  which 
employs  a  finite-size  pencil  beam  (FSPB)  algorithm,  agreed  with  the  Monte  Carlo  simulations 
and  measurements  to  within  4%  in  a  cylindrical  water  phantom  with  various  hypothetical  target 
shapes.  Discrepancies  of  more  than  5%  (relative  to  the  prescribed  target  dose)  in  the  target  region 
and  over  20%  in  the  critical  structures  were  found  in  some  IMRT  patient  calculations.  The  FSPB 
algorithm  as  implemented  in  the  Corvus  system  is  adequate  for  homogeneous  phantoms  (such  as 
prostate)  but  may  result  in  significant  under-  or  over-estimation  of  the  dose  in  some  cases  involving 
heterogeneities  such  as  the  air-tissue,  lung-tissue  and  tissue-bone  interfaces. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version;  see  www.iop.org) 


1,  Introduction 

For  the  last  few  years ,  extensive  research  has  been  carried  out  to  develop  conformal  radiotherapy 
using  computer-controlled  linear  accelerators  equipped  with  multileaf  collimators  (MLC) 
(Boesecke  et  al  1988,  Leibel  et  al  1992,  LoSasso  et  al  1993,  Powlis  et  al  1993,  Mageras 
et  al  1994,  Brewster  et  al  1995,  Fraass  et  al  1995,  McShan  et  al  1995,  Yu  et  al  1995).  More 
recently,  intensity  modulated  radiotherapy  (IMRT)  has  been  developed  (Brahme  1988,  Convery 
and  Rosenbloom  1992,  Webb  1992,  Boyer  et  al  1997)  and  implemented  (Ling  et  al  1996, 
Boyer  et  al  1 998)  that  uses  computer-controlled  modulation  of  x-ray  fields  by  the  MLC.  It  is 
anticipated  that  conformal  radiotherapy  will  provide  radiation  oncologists  with  a  significantly 
improved  tool  to  deliver  high  doses  of  ionizing  radiation  to  some  tumour  sites  while  reducing 
doses  to  adjacent  normal  tissue  below  levels  to  which  they  are  unavoidably  exposed  by  currently 
available  techniques.  Thus,  acute  and  chronic  toxicity  associated  with  treatment  of  a  tumour 
volume  by  radiation  may  be  significantly  reduced  or  delayed  for  certain  sites  of  malignant 
presentations. 
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The  use  of  conformal  radiotherapy,  especially  with  the  IMRT  technique,  is  a  major 
departure  from  the  way  radiotherapy  is  currently  delivered.  Although  the  use  of  MLCs 
provides  the  possibility  of  achieving  better  dose  distributions  conformed  to  tumour  targets, 
it  also  increases  the  complexity  of  treatment.  The  sequences  of  leaf  movement  and  their 
associated  effects  on  the  dose  delivered  to  the  patient  may  vary  significantly  depending  on  the 
accelerator  and  the  MLC  design.  Important  factors  include  the  variation  of  the  accelerator 
head  scatter  component  in  the  MLC^colli mated  beam  (Convery  and  Webb  1997),  the  amount 
of  photon  leakage  through  the  leaves  (Wang  et  al  1996,  Webb  1997,  Holmes  et  al  1997),  the 
scatter  from  the  leaf  ends,  the  ‘tongue  and  groove’  effect  (Chui  et  al  1994,  Wang  etal  1996), 
the  effect  of  back-scattered  photons  from  the  moving  jaws  and  MLC  leaves  on  the  monitor 
chamber  signal  (Hounsell  1998).  Traditionally,  patient  dose  calculations  in  radiotherapy  have 
been  based  on  correcting  measured  dose  distributions.  New  dose  calculation  algorithms  have 
been  developed  to  predict  the  patient  dose  from  ‘first  principles’  using  a  model  of  radiation 
transport  (Mackie  etal  1995).  Comparisons  of  the  traditional  photon  algorithms  and  the  newer 
ones  have  been  reviewed  by  Wong  and  Purdy  (1990),  Cunningham  and  Battista  (1995)  and 
Mackie  et  al  (1996).  Due  to  the  lack  of  electron  transport,  the  conventional  dose  calculation 
algorithms  often  failed  to  predict  the  dose  distribution  accurately  near  inhomogeneities  (Mackie 
etal  1996,  Mohan  1997,  DeMarco  1998,  Wang  era/  1998,  Mae/ a/  1999).  Furthermore, 
the  inverse-planning  algorithms  for  beam  optimization  have  all  used  approximations  to  speed 
up  the  dose  computation  that  may  introduce  significant  uncertainty  in  the  calculated  dose 
distributions,  especially  in  the  presence  of  heterogeneities.  When  simple  source  models  are 
used  in  the  dose  computation,  the  correlation  between  the  calibrated  reference  dose  and  the 
dose  related  to  a  beam  segment  may  be  lost.  All  the  above  imply  a  potential  problem  with  the 
prediction  of  the  dose  distributions  in  a  patient  for  an  IMRT  treatment. 

Oldham  and  Webb  (1997)  reported  differences  in  excess  of  10%  in  the  absolute  dose 
between  the  optimization  dose  calculations  and  measurements  (using  film)  of  fields  delivered 
by  a  dynamic  MLC.  The  differences  were  attributed  partially  to  the  nonlinearity  of  dose  per 
monitor  unit  (MU)  for  small  MU  deliveries  (the  actual  dose  delivered  per  MU  increased  by 
more  than  1 0%  from  MU  =  20  to  MU  =  1 ).  Wang  etal  {\ 996)  reported  that  for  the  Memorial 
Hospital’s  dynamic  MLC  delivery  process,  the  discrepancies  between  the  calculated  dose  and 
the  measured  dose  were  in  excess  of  5%  if  various  effects  related  to  the  MLC  construction, 
such  as  accelerator  head  scatter,  were  not  properly  accounted  for.  The  uncertainty  in  the 
doses  calculated  by  a  conventional  dose  calculation  algorithm  was  5-10%  in  the  presence  of 
heterogeneities  (Mohan  1997).  Our  recent  Monte  Carlo  results  were  consistent  with  these 
findings  (Ma  et  al  1999). 

The  purpose  of  this  work  was  to  verify  the  accuracy  of  the  IMRT  dose  distributions  from 
a  commercial  treatment  planning  optimization  system  (Corvus,  Nomos  Corp.,  Sewickley,  PA) 
using  Monte  Carlo  simulations.  We  have  used  the  EGS4/BEAM  Monte  Carlo  code  system 
(Nelson  et  al  1985,  Rogers  et  al  1995a,  b)  to  simulate  the  clinical  photon  beams  from  two 
linear  accelerators,  Varian  Clinac  2100C  and  2300C/D  (Varian  Oncology,  Palo  Alto,  CA).  The 
EGS4/DOSXYZ  code  (Rogers  et  al  1995a,  Ma  et  al  1995)  was  modified  to  compute  the  dose 
in  a  patient  (a  three-dimensional  phantom  built  from  the  CT  data)  for  multiple  fixed-gantry 
fields.  Phantom  measurements  were  performed  to  commission  the  Monte  Carlo  system.  The 
dose  distributions  in  the  experimental  phantoms  and  in  the  patients  were  calculated  and  used 
to  verify  the  optimized  treatment  plans  generated  by  the  Corvus  inverse-planning  system.  In 
the  following  sections,  we  will  describe  the  dose  calculations  in  the  inverse-planning  system 
and  the  Monte  Carlo  simulations.  We  will  show  the  dose  distributions  for  homogeneous  and 
heterogeneous  phantoms  and  discuss  the  effect  of  material  density  and  atomic  number  on  the 
final  dose  calculations. 
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2.  Materials  and  method 

2. 1.  Treatment  planning  and  dose  measurement 

The  treatment  planning  optimization  system  used  in  this  work  is  the  Corvus  inverse-planning 
system.  The  inverse  planning  process  employs  a  simulated  annealing  optimization  algorithm 
(Webb  1992).  The  dose  calculation  algorithm  used  for  the  inverse  planning  process  is  a  finite- 
size  pencil  beam  (FSPB)  algorithm,  which  uses  predetermined  beamlet  dose  distributions.  The 
beamlet  distributions  were  derived  from  measured  dose  distributions  and  normalized  to  produce 
consistent  output  factors  for  various  field  sizes.  The  patient  inhomogeneity  correction  is  made 
by  ‘stretching’  the  beamlet  distribution  proportionally  based  on  the  equivalent  pathlength. 
Some  details  of  the  dose  calculations  have  been  described  by  Holmes  et  al  (1998)  and  Boyer 
et  al  (1999).  The  monitor  unit  calculations  for  a  ‘step-and-shoot’  leaf  sequence  algorithm 
have  been  discussed  by  Boyer  et  al  (1999).  The  effect  of  the  leaf  leakage  is  accounted  for 
in  Corvus  by  reducing  the  original  beamlet  weight  by  the  same  amount  as  the  accumulated 
leakage  for  a  given  leaf  sequence.  This  method  works  well  for  those  beamlets  whose  weights 
are  greater  than  the  estimated  leakage.  If  the  beamlet  weights  are  smaller  than  the  leakage 
this  method  will  underestimate  the  dose  as  any  remaining  leakage  effect  cannot  be  accounted 
for  and  the  initial  dose  calculation  does  not  include  any  leaf  leakage  effect  (which  is  unknown 
before  a  leaf  sequence  is  generated  for  the  field).  The  Corvus  system  has  been  commissioned 
for  clinical  IMRT  treatment  planning  (Xing  et  al  1999)  and  used  for  treating  head  and  neck 
patients  with  IMRT  (Boyer  et  al  1998). 

In  order  to  test  the  Corvus  system,  inverse  plans  were  computed  for  various  target  shapes 
placed  in  the  centre  of  a  cylindrical  water  phantom  having  a  diameter  of  30  cm.  Comparisons 
of  the  Corvus  dose  distributions  with  measurements  have  been  reported  in  detail  by  Boyer  et  al 
(1999).  In  this  work,  we  have  computed  the  dose  distributions  for  these  plans  using  Monte 
Carlo  methods  (see  descriptions  below).  For  completeness,  we  briefly  describe  the  plans  and 
the  measurements  below.  Inverse  plans  were  computed  for  different  hypothetical  targets  with 
different  numbers  of  beams  directed  toward  the  axis  of  the  cylinder  at  the  centre  of  the  target  and 
spaced  at  equal  angles.  The  treatments  were  delivered  using  a  dynamic  multileaf  collimator 
(Varian  Oncology  Systems,  Palo  Alto,  CA).  The  leaf  sequences  computed  within  the  Nomos 
software  were  written  into  files  with  formats  conformed  to  the  requirements  of  the  Varian 
digital  control  software.  The  leaf  sequences  were  delivered  using  the  monitor  units  calculated 
by  the  Corvus  system.  The  absolute  dose  delivered  by  the  leaf  sequences  was  measured  in  the 
30  cm  diameter  cylindrical  water  phantom  using  a  0.147  cm^  ionization  chamber  (Wellhofer 
Dosimetric,  Schwartzenbruck,  Germany)  following  the  AAPM  TG-21  protocol  (AAPM  1983). 
No  corrections  were  made  for  the  variation  in  the  chamber  displacement  effect,  which  depends 
on  the  dose  gradient  at  the  measurement  point  and  the  chamber  diameter.  This  may  introduce 
an  up  to  2%  uncertainty  in  the  measured  dose  for  the  6  mm  diameter  chamber  used  (the  dose 
gradient  was  5-8%  per  centimetre  at  some  measurement  points).  The  chamber  positioning 
uncertainty  was  about  0.1  cm.  The  overall  uncertainty  in  the  measured  dose  was  estimated  to 
be  about  3%  (lor). 

2.2.  The  Monte  Carlo  simulation 

We  have  used  BEAM  and  DOSXYZ  (Rogers  etal  1995a,  b,  Ma  etal  1995)  Monte  Carlo  codes 
for  the  accelerator  head  simulation  and  dose  calculation  in  the  patient  respectively.  Both  codes 
were  EGS4  (Electron  Gamma  Shower  version  4,  Nelson  etal  (1985))  user  codes,  running  under 
the  UNIX  operating  system,  developed  through  the  Omega  project  for  Monte  Carlo  treatment 
planning  dose  calculations  (Mackie  et  al  1994).  Detailed  descriptions  of  the  software  can  be 
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found  from  Rogers  et  al  (1995a,  b).  A  detailed  description  of  the  clinical  implementation  of 
the  Monte  Carlo  code  system  was  given  in  a  previous  publication  (Ma  et  al  1999). 

Two  types  of  clinical  linear  accelerators  were  simulated  for  the  clinical  implementation 
of  Monte  Carlo  dose  verification  in  this  work:  Varian  Clinac  2100C  and  2300C/D  (Varian 
Oncology  Systems,  Palo  Alto,  CA).  The  dimensions  and  materials  for  the  accelerator 
components  were  incorporated  according  to  the  manufacturer’s  specifications.  Electron 
beams  emerging  from  the  vacuum  exit  window  were  assumed  to  be  monoenergetic  and 
monodirectional  with  a  beam  radius  of  0. 1-0.2  cm.  These  were  found  to  be  reasonable 
assumptions  to  achieve  an  acceptable  dose  calculation  accuracy  of  about  2%  of  the  dose 
maximum  (Anax)  anywhere  in  the  phantom  for  clinical  radiotherapy  applications  (Kapur  etal 
1998,  Ma  1998,  Ma  et  al  1999).  We  obtained  accurate  phase-space  data  for  photon  beams 
with  nominal  energies  of  4,  6  and  1 5  MV.  The  energy  cut-offs  for  electron  transport  in  the 
accelerator  simulation  (ECUT  and  AE)  700  keV  (kinetic  +  rest  mass)  and  for  photon  transport 
(PCUT  and  AP)  10  keV.  The  maximum  fractional  energy  loss  per  electron  step  (ESTEPE)  was 
0.04.  The  bremsstrahlung  splitting  and  Russian  roulette  options  were  implemented  for  photon 
beam  simulations  (Rogers  etal  1995b).  The  ICRU  recommended  compositions  and  stopping 
power  values  were  used  for  the  materials  used  in  the  accelerator  simulations  (ICRU  1984). 
The  phase-space  data  were  scored  at  a  plane  immediately  above  the  photon  jaws.  The  number 
of  particles  in  a  photon  beam  file  was  about  50  million.  Field  shaping  by  photon  jaws,  blocks 
and  the  MLC  was  further  simulated  using  BEAM  and  the  phase-space  data  could  be  stored 
temporarily  or  used  directly  for  dose  calculations. 

The  DOSXYZ  code  was  designed  for  dose  calculations  in  a  3D  rectilinear  voxel  (volume 
element)  geometry  (Ma  et  al  1995).  Voxel  dimensions  were  completely  variable  in  all  three 
directions.  Every  voxel  could  be  assigned  to  a  different  material.  The  cross-section  data 
for  the  materials  used  were  available  in  a  pre-processed  PEGS4  cross-section  data  file.  The 
mass  density  of  the  material  in  a  DOSXYZ  calculation  was  varied  based  on  the  patient’s 
CT  data,  although  the  density  effect  corrections  for  the  stopping  powers  of  the  material 
remain  unchanged  (Ma  et  al  1999).  The  voxel  dimensions  and  materials  were  defined  in 
a  DOSXYZ  input  file  together  with  the  transport  parameters  such  as  ECUT,  PCUT,  ESTEPE 
and  the  parameters  required  by  PRESTA  (Bielajew  and  Rogers  1988).  The  phase-space 
data  obtained  from  a  BEAM  simulation  was  used  as  a  source  input  with  variable  source 
positions  and  beam  incidence  angles.  Dose  contributions  from  different  beam  components 
were  selectively  calculated  based  on  the  particle  charge  or  the  LATCH  settings  specified  in 
the  BEAM  simulation.  DOSXYZ  produced  a  data  file  that  contained  geometry  specifications 
such  as  the  number  of  voxels  in  all  three  directions  and  their  boundaries  as  well  as  the  dose 
values  and  the  associated  (la)  statistical  uncertainties  in  the  individual  voxels. 

The  DOSXYZ  code  has  been  modified  to  read  the  MLC  leaf  sequence  files  for  IMRT 
treatment.  For  this  study,  the  Varian  dynamic  MLC  leaf  sequence  files  (G-version)  generated 
using  a  step-and-shoot  leaf  sequencing  algorithm,  were  used.  The  monitor  units  for  each 
leaf  sequence  were  integrated  into  a  two-dimensional  map  (grids).  The  dimensions  of  the 
map  at  the  isocentre  (100  cm  from  the  virtual  point  source  position)  are  40  cm  in  both  the  x 
and  y  directions  and  the  grid  sizes  were  0.1  cm.  The  monitor  units  for  the  open  areas  were 
accumulated  for  each  of  the  pixels,  while  for  the  closed  areas  a  fraction  of  the  monitor  units 
were  accumulated  based  on  the  measured  MLC  leakage  factor  for  each  beam  energy.  Thus, 
the  integrated  intensity  map  has  included  the  averaged  leaf  leakage  effect  but  ignored  the 
influence  of  the  leaf  shape  (tongue  and  groove)  and  the  variation  of  the  spatial  and  spectral 
distributions  due  to  photon  attenuation  and  scattering  in  the  MLC  leaves.  This  is  considered  to 
be  a  reasonable  approach  since,  unless  the  effect  of  organ  motion  and  patient  set-up  uncertainty 
has  been  accounted  for  properly,  the  averaged  leaf  leakage  is  a  more  realistic  quantity  than 
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the  explicit  leakage  through  the  leaf  and  the  tongue  and  groove  for  estimating  the  overall  dose 
in  a  patient  during  the  whole  treatment  course.  The  transmission  factor  for  the  areas  under 
the  photon  jaws  was  assumed  to  be  zero.  The  dimensions  of  the  photon  jaw  opening  were  set 
according  to  the  actual  treatment  set-ups.  Further  modifications  were  made  to  simulate  several 
gantry  angles  in  the  same  run.  The  leaf  sequence  files  were  read  sequentially  and  the  gantry 
angle  was  changed  automatically  after  the  simulation  of  a  photon  field  was  completed. 

During  the  Monte  Carlo  simulation,  the  weight  of  a  phase-space  particle  was  altered  based 
on  the  value  of  the  pixel  in  the  intensity  map  through  which  the  particle  was  travelling.  To 
improve  the  simulation  efficiency,  two  variance  reduction  techniques  were  implemented  in 
the  DOSXYZ  code,  particle  splitting  and  Russian  roulette  (Rogers  and  Bielajew  1990).  The 
former  was  applied  to  the  particles  with  weight  greater  than  unity  and  the  latter  was  applied 
to  the  particles  with  weight  smaller  than  unity.  The  splitting  and  Russian  roulette  processes 
were  implemented  in  such  a  way  that  the  resulting  particles  would  have  identical  weighting 
factors.  Using  uniform  weighting  factors  will  generally  improve  the  statistical  uncertainty  of 
a  simulation. 

For  patient  dose  calculations,  the  simulation  phantom  was  built  from  the  patient’s  CT  data 
with  up  to  128  X  128  X  128  voxels  (uniform  in  any  dimensions).  The  side  of  a  voxel  varied 
from  0.2  to  0.4  cm.  The  EGS4  transport  parameters  were  set  to  ECUT  =  AE  =  700  keV, 
PCUT  =  AP  =  10  keV  and  ESTEPE  =  0.04.  The  number  of  particle  histories  simulated 
ranged  from  300  million  to  1  billion  for  up  to  nine  gantry  angles  for  an  IMRT  treatment.  The 
lor  statistical  uncertainty  in  the  dose  was  generally  1-2%  of  the  Dmax  value.  The  CPU  time 
required  for  an  IMRT  simulation  was  3-15  h  on  a  single  Pentium  III  450  MHz  PC  depending 
on  the  beam  energy  and  field  sizes. 


3.  Results 

3.L  Comparisons  in  homogeneous  water  phantoms 

Both  the  inverse-planning  system  and  the  Monte  Carlo  dose  calculation  system  have  been 
commissioned  for  routine  clinical  applications  (Xing  et  al  1999,  Ma  1998,  Ma  et  al  1999). 
Both  systems  regenerated  the  PDD  curves  and  the  dose  profiles  at  various  depths  to  within 
2%  of  the  dose  maximum  values  of  the  measured  beam  data  for  various  field  sizes  and  source 
to  surface  distances  (SSD).  The  dose  distributions  given  by  the  Corvus  system  agreed  with 
measurements  to  about  3%  in  a  cylindrical  water  phantom  with  various  hypothetical  target 
shapes  (Boyer  et  al  1999).  Figure  1  shows  the  dose  distributions  in  a  30  cm  diameter  water 
cylinder  irradiated  by  seven  6  MV  intensity-modulated  photon  beams  calculated  by  Corvus 
and  by  the  Monte  Carlo  simulations.  The  critical  structure  is  in  the  centre  of  the  water  cylinder 
immediately  next  to  the  C-shaped  target.  A  beamlet  size  of  1  cm  x  1  cm  was  used  in  the 
optimization  calculation.  The  leaf  sequencing  algorithm  used  20  intensity  levels.  The  isodose 
lines  calculated  by  both  systems  agreed  within  about  4%  or  a  shift  in  isodose  lines  within  about 
0.3  cm.  These  small  discrepancies  were  thought  to  be  partially  due  to  the  approximations  used 
by  Corvus  in  accounting  for  the  effect  of  MLC  leaf  leakage.  Considering  that  the  average  leaf 
leakage  is  about  1 .5%  for  a  6  MV  beam  and  the  total  MU  required  to  deliver  an  IMRT  field  is 
about  three  times  more  than  that  required  by  a  conventional  field,  the  estimated  accumulated 
leaf  leakage  beneath  a  completely  closed  leaf  pair  is  about  4.5%.  Thus,  a  few  per  cent  of  the 
observed  dose  difference  may  be  attributed  to  the  difference  in  leaf  leakage  implementation 
between  Corvus  and  Monte  Carlo. 
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Figure  1.  Comparison  of  6  MV  photon  beam  (seven  co-planar  fields  at  26,  77,  180,  231,  283  and 
334°  gantry  angles)  dose  distributions  in  a  water  cylinder  calculated  by  Monte  Carlo  (a)  and  Corvus 
{b).  The  differences  between  the  Corvus  calculations  and  measurements  using  an  ion  chamber  are 
shown  for  certain  locations  (arrows). 


3.2.  Comparisons  in  CT phantoms 

In  this  work,  we  have  compared  the  dose  distributions  calculated  by  Monte  Carlo  and  the 
Corvus  FSPB  algorithm  for  various  treatment  sites  to  identify  potential  treatment  situations 
that  may  benefit  from  dose  accuracy  improvements.  In  the  following  examples,  we  show  two 
typical  IMRT  treatment  plans  computed  by  the  Corvus  inverse  planning  system  and  verified 
by  the  Monte  Carlo  system. 

3.2.1.  Prostate.  To  explore  the  effect  of  photon  and  electron  transport  on  IMRT  dose 
calculations,  we  show  in  figure  2  an  IMRT  prostate  treatment  plan  calculated  by  Monte  Carlo 
simulation  and  the  Corvus  system.  The  plan  was  generated  using  the  Corvus  system  for  15  MV 
photon  beams  with  nine  gantry  angles  (20,  60,  100,  140,  220,  260,  300  and  340°).  The  beam 
intensity  was  modulated  using  a  Varian  dynamic  MLC  with  80  leaves.  In  both  calculations, 
the  isodose  lines  represent  the  absolute  dose  values  in  the  patient.  It  can  be  seen  that  the  dose 
values  in  the  target  (the  prostate)  agreed  very  well  between  the  Corvus  calculations  and  the 
Monte  Carlo  simulations.  Similar  results  were  found  for  six  other  prostate  cases  compared 
(not  shown).  The  difference  in  the  average  dose  to  the  target  volume  was  1.6%  between 
Monte  Carlo  and  Corvus  and  the  maximum  dose  difference  in  the  target  dose  was  3.4%. 
This  confirms  that  both  calculation  algorithms  can  predict  dose  distributions  in  homogeneous 
phantoms  accurately.  The  dose  values  in  the  nearby  critical  structures  also  agreed  to  3-7% 
of  the  prescribed  target  dose.  For  the  case  shown  in  figure  2(c),  the  maximum  difference  was 
about  3  Gy  in  the  rectum.  However,  other  cases  also  showed  the  maximum  differences  in  the 
bladder.  These  differences  may  be  considered  to  be  clinically  acceptable. 

The  dose  values  in  the  regions  near  the  bony  structures  sometimes  showed  a  difference  of 
a  few  per  cent  between  the  Corvus  and  the  Monte  Carlo  calculations.  These  discrepancies  may 
be  partially  explained  by  the  effect  of  electron  backscatter  from  the  bone  (high  atomic  number 
and  high  density)  to  the  soft  tissue,  which  was  accurately  accounted  for  by  the  Monte  Carlo 


MC  verification  of  IMRT  dose  distributions 


2489 


Figure  2.  Comparison  of  15  MV  photon  beam  (nine  co-planar  fields)  dose  distributions  for  a 
prostate  treatment  calculated  by  Monte  Carlo  (a)  and  Corvus  (h).  The  isodose  lines  are  112, 70.0, 
56.1,48.9,35.0, 27.8,21.1  and  13.9  Gy  respectively  in  each  figure,  (c)  The  dose  volume  histograms 
as  calculated  by  Monte  Carlo  (full  curves)  and  Corvus  (broken  curves). 


simulation  but  not  in  a  FSPB  algorithm.  There  have  been  several  Monte  Carlo  studies  in  the 
literature  showing  similar  discrepancies  between  the  Monte  Carlo  algorithm  and  the  correction 
based  calculation  algorithms  in  dose  build-up  or  build-down  regions  near  air  cavities,  lung  and 
large  bony  structures  (DeMarco  et  al  1998,  Mohan  1997,  Wang  et  al  1996,  Ma  et  al  1999).  It 
should  be  mentioned  that  the  difference  in  the  dose  to  the  bone  was  partially  due  to  the  fact 
that  the  conventional  dose  calculation  algorithms  usually  used  water  as  the  phantom  material 
and  the  inhomogeneity  corrections  were  computed  using  varying  electron  density  (based  on 
the  CT  numbers)  while  the  Monte  Carlo  algorithm  used  different  materials  such  as  air,  tissue, 
lung  and  bone  with  varying  mass  density  calculated  from  the  CT  data.  If  we  convert  the  dose 
to  the  bone  material  to  the  dose  to  tissue  using  the  stopping  power  ratio  for  bone  to  tissue 
(assuming  the  same  electron  energy  fluence)  the  dose  in  the  bone  regions  will  be  about  3.5% 
higher  for  soft  bone  and  about  10%  higher  for  compact  bone  (ICRU  1984,  Siebers  etal  2000). 
We  will  discuss  this  issue  further  in  section  4.  Anyway,  the  dose  in  the  surrounding  tissue 
regions  should  not  be  affected  by  the  conversion,  which  was  mainly  caused  by  the  perturbation 
of  the  electron  fluence  by  the  nearby  inhomogeneous  anatomy. 

5.2.2.  Vertebra,  Figure  3  shows  the  dose  distributions  for  the  treatment  of  the  vertebra 
calculated  by  both  Monte  Carlo  simulation  and  the  Corvus  system.  The  plan  was  generated 
using  the  Corvus  system  for  15  MV  photon  beams  with  nine  co-planar  gantry  angles  (20,  55, 
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Figure  3.  Dose  distributions  for  the  treatment  of  the  vertebra  calculated  by  Monte  Carlo  (a)  and 
by  Corvus  (h)  for  15  MV  photons  (nine  co-planar  fields).  The  isodose  lines  are  17.6,  15.6,  13.7, 
1 1.7,  9.8,  7.8,  5.9,  3.9  and  2.0  Gy  respectively  in  each  figure,  (c)  The  dose  volume  histograms  as 
calculated  by  Monte  Carlo  (full  curves)  and  Corvus  (broken  curves)  for  the  target  and  the  spinal 
cord. 


90,  140,  180,  220,  260,  300  and  340°).  The  intensity  was  modulated  using  a  Varian  dynamic 
MLC  with  80  leaves.  The  prescribed  target  dose  was  18  Gy.  The  maximum  dose  in  the  target 
showed  good  agreement  between  Corvus  and  Monte  Carlo  (figure  3(c)).  The  Monte  Carlo 
dose  distribution  showed  slightly  better  target  coverage  than  the  Corvus  dose  distribution  (a 
2  Gy  difference  in  the  minimum  target  dose).  Because  the  cord  was  sometimes  immediately 
next  to  the  target  region  the  maximum  cord  dose  was  expected  to  be  equal  to  or  higher  than  the 
minimum  target  dose.  This  was  confirmed  by  the  Monte  Carlo  simulations  (see  figure  3(c)). 
In  the  regions  near  large  bony  structures  (such  as  the  cord)  differences  of  more  than  20%  of  the 
prescribed  target  dose  could  be  seen  between  the  Corvus  calculation  (10  Gy)  and  the  Monte 
Carlo  simulation  (14  Gy).  The  difference  in  the  dose  to  the  cord  was  thought  to  be  due  in 
part  to  electron  scattering  from  the  surrounding  bone,  which  could  not  be  modelled  properly 
using  the  FSPB  algorithm.  Another  possible  reason  might  be  due  to  the  implementation  of  the 
heterogeneity  and  leaf  leakage  corrections  in  the  FSPB  model.  Although  the  photon  beams 
were  optimized  to  avoid  the  cord,  electrons  could  reach  the  cord  and  the  dose  to  the  cord  could 
be  enhanced  due  to  the  high-density  material  surrounding  it  and/or  photon  leaf  leakage  that  is 
not  included  in  the  dose  calculation  during  the  inverse  planning  process.  Further  studies  are 
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Figure  4.  Dose  distributions  for  different  tissue  types  and  material  densities:  (a)  tissue  and  bone 
with  variable  density  (thick  line)  and  tissue  with  unity  density  (thin  line);  {b)  tissue  and  bone  with 
variable  density  (thick  line)  and  tissue  with  variable  density  (thin  line);  and  (c)  tissue  and  bone  with 
variable  density  (thick  line)  and  tissue  with  unity  density  and  bone  with  10  g  cm“^  density  (thin 
line).  The  phantom  geometry  and  beam  arrangements  are  the  same  as  in  figure  3.  The  isodose 
lines  are  given  as  10,  30,  50,  70  and  80%  of  the  prescribed  target  dose. 


needed  to  understand  these  differences  if  more  access  to  the  FSPB  and  leaf  sequence  algorithms 
in  the  Corvus  system  is  available. 

4.  Discussion 

Several  important  factors  may  affect  the  Monte  Carlo  calculated  dose  distributions  and  the 
DVH  curves.  First,  the  isodose  lines  and  the  DVH  curves  are  affected  by  the  materials  used  in 
the  patient  CT  phantom,  i.e.  whether  we  plot  dose  to  tissue  only  or  dose  to  any  material  (such 
as  air,  tissue  or  bone).  It  seems  reasonable  that  our  previous  experience  was  based  on  dose  to 
tissue  (or  dose  to  water,  the  difference  between  the  two  is  within  1%)  and  therefore  the  dose 
values  should  be  expressed  as  dose  to  tissue.  However,  it  can  also  be  argued  that  the  real  dose 
to  the  biological  material  such  as  bone  should  be  given  whenever  possible.  Only  in  this  way 
can  the  relationship  between  the  ‘old’  practice  and  new  experience  be  established. 

To  understand  the  effect  of  the  conversion  of  the  dose  to  different  materials,  we  show 
in  figure  4  the  dose  distributions  calculated  using  Monte  Carlo  with  different  materials  and 
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Figure  5.  Dose  volume  histograms  for  the  GTV  in  figure  3.  Full  circles  are  obtained  using  a 
CT  phantom  consisting  of  air,  tissue  and  bone  of  variable  density  with  dose  calculated  for  each 
material.  Open  circles  are  obtained  using  the  same  phantom  with  dose  converted  to  tissue  using  the 
stopping  power  ratios  of  tissue  to  bone  or  tissue  to  air.  Open  squares  are  obtained  using  a  phantom 
consisting  of  only  tissue  of  variable  density. 


density  configurations  for  the  IMRT  case  shown  in  figure  3.  In  figure  4(«),  the  dose  distribution 
was  calculated  using  tissue  with  unit  density  and  air.  This  phantom  should  not  show  any 
heterogeneity  effect  due  to  the  change  in  tissue  (or  bone)  densities.  In  figure  A{b),  the  dose 
distribution  was  calculated  using  air  and  tissue  with  variable  density  converted  from  the  CT 
data.  This  is  similar  to  the  phantom  used  by  the  Corvus  system.  In  figure  4(c),  the  dose 
distribution  was  calculated  with  air,  unit  density  tissue  and  10  g  cm“^  density  bone.  For 
comparison,  the  dose  distribution  calculated  using  air,  tissue  and  bone  with  proper  densities 
converted  from  the  CT  data  is  shown  in  each  figure.  The  isodose  curves  were  computed  by 
normalizing  the  dose  values  to  the  prescribed  target  dose.  The  differences  in  the  isodose  lines 
in  figure  4(a)  and  {b)  are  small  (2-3%).  This  is  also  clearly  shown  in  figure  5  for  the  target 
DVHs.  The  difference  between  dose  to  tissue  and  dose  to  bone  would  be  about  3.5%  for  soft 
bone  and  about  10%  for  hard  bone  (Siebers  et  al  2000)  if  converted  using  the  stopping  power 
ratios  for  tissue  to  bone  (or  soft  bone)  at  these  beam  qualities,  assuming  the  same  electron 
energy  fluence.  Clearly,  such  conversion  is  not  equivalent  to  performing  the  Monte  Carlo 
simulation  using  a  water  (or  tissue)  phantom  with  variable  mass  (or  electron)  density.  It  is 
known  that  electron  backscattering  from  the  high  atomic  number  materials  may  perturb  the 
dose  in  tissue  near  the  tissue-bone  interface.  This  effect  is  less  significant  when  dose  values  are 
averaged  over  course  scoring  volumes  (0.3-0.4  cm  voxels).  In  figure  4(c),  the  density  of  bone 
was  artificially  increased  to  10  g  cm“^  which  caused  significant  attenuation  of  the  beams  and 
therefore  altered  the  doses  behind  the  bones.  The  effects  on  the  surface  doses  are  smaller  than 
for  the  high  dose  regions.  The  maximum  dose  for  the  artificially  high-density  bone  geometry 
is  about  15%  lower  than  that  for  the  phantom  with  normal  material  density. 
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Although  Monte  Carlo  dose  calculations  are  time-consuming  it  is  possible  to  use  Monte 
Carlo  calculated  dose  distributions  to  verify  the  IMRT  treatment  plans  in  order  to  detect  any 
cases  where  the  FSPB  dose  calculation  algorithm  as  implemented  in  Corvus  may  fail  to  predict 
the  dose  perturbation  effect  near  inhomogeneities.  Monte  Carlo  simulation  may  also  be  used 
directly  for  IMRT  beamlet  distribution  calculation  as  a  practical  solution  to  this  problem.  Our 
experiences  show  that  a  factor  of  two  to  three  more  Monte  Carlo  particle  histories  are  needed 
for  an  IMRT  treatment  simulation  compared  with  a  conventional  photon  treatment  simulation 
to  achieve  the  same  statistical  uncertainty.  This  is  because  more  monitor  units  are  needed  to 
deliver  intensity  modulated  photon  fields;  more  particles  will  be  simulated  in  a  Monte  Carlo 
calculation  but  many  of  them  will  be  stopped  by  the  MLC  leaves.  Therefore,  the  CPU  time 
per  photon  history  for  an  IMRT  simulation  is  less  than  that  for  a  conventional  field.  Using  the 
existing  computing  power  (8  CPUs)  of  the  Corvus  system,  the  calculation  time  for  a  typical 
‘inverse  plan’  would  be  increased  from  the  current  0.5-1  h  to  2-4  h  with  two  Monte  Carlo 
calculations.  The  pre-optimization  dose  calculation  will  provide  the  beamlet  distributions  for 
the  optimization  process,  which  take  into  account  the  effect  of  the  accelerator  head  scatter  and 
inhomogeneous  anatomy  of  the  patient.  The  post-optimization  dose  calculation  will  include  the 
effects  due  to  leaf  leakage,  leaf  scatter  and  photon  backscatter  into  the  monitor  chamber  after 
the  sequence  of  MLC  leaf  movement  and  jaw  positions  has  been  generated.  Further  studies  are 
under  way  on  a  Monte  Carlo  dose  calculation  based  inverse  planning  system  (Pawlicki  et  al 
1999,  Ma  etal  2000). 

5.  Summary 

We  have  implemented  a  Monte  Carlo  system  for  routine  radiotherapy  treatment  planning 
dose  calculations.  In  our  previous  publications  we  have  shown  that  Monte  Carlo  simulations 
agreed  with  measurements  to  within  2%  for  various  clinical  beam  set-ups  in  homogeneous 
and  heterogeneous  phantoms.  Based  on  these  results,  we  have  moved  one  step  nearer  to  using 
the  Monte  Carlo  simulations  to  verify  the  IMRT  dose  distributions  computed  by  the  Corvus 
system,  which  employs  a  FSPB  algorithm  for  beamlet  dose  calculations,  assuming  that  the 
Monte  Carlo  simulations  are  accurate  at  the  2%  level  for  patient  phantoms  built  from  CT  data. 

Our  results  showed  that  the  FSPB  algorithm  was  adequate  for  most  of  the  IMRT  cases 
where  the  target  was  not  immediately  adjacent  to  the  critical  structures.  However,  the  FSPB 
algorithm  may  not  accurately  predict  the  dose  distributions  in  and  near  inhomogeneities  in 
some  cases.  The  dose  in  the  target  volume  calculated  by  the  Corvus  system  differed  from  the 
Monte  Carlo  results  by  more  than  5%,  while  the  dose  to  the  critical  organ  differed  by  more 
than  20%  of  the  prescribed  target  dose  for  a  few  cases.  This  suggests  that,  for  such  cases, 
more  accurate  dose  calculation  algorithms  than  that  currently  implemented  in  Corvus  should 
be  used  for  intensity-modulated  radiotherapy  treatment  planning. 
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Abstract  Monte  Carlo  simulations  were  employed  to  study  the  characteristics  of  the  electron 
beams  of  a  clinical  linear  accelerator  in  the  presence  of  1.5  and  3.0  T  transverse  magnetic  fields  and 
to  assess  the  possibility  of  using  magnetic  fields  in  conjunction  with  modulated  electron  radiation 
therapy  (MERT).  The  starting  depth  of  the  magnetic  field  was  varied  over  several  centimetres.  It 
was  found  that  peak  doses  of  as  much  as  2.7  times  the  surface  dose  could  be  achieved  with  a  1.5  T 
magnetic  field.  The  magnetic  field  was  shown  to  reduce  the  80%  and  20%  dose  drop-off  distance 
by  50%  to  80%.  The  distance  between  the  80%  dose  levels  of  the  pseudo-Bragg  peak  induced 
by  the  magnetic  field  was  found  to  be  extremely  narrow,  generally  less  than  1  cm.  However, 
by  modulating  the  energy  and  intensity  of  the  electron  fields  while  simultaneously  moving  the 
magnetic  field,  a  homogeneous  dose  distribution  with  low  surface  dose  and  a  sharp  dose  fall-off 
was  generated.  Heterogeneities  are  shown  to  change  the  effective  range  of  the  electron  beams, 
but  not  eliminate  the  advantages  of  a  sharp  depth-dose  drop-off  or  high  peak-to-surface  dose  ratio. 
This  suggests  the  applicability  of  MERT  with  magnetic  fields  in  heterogeneous  media.  The  results 
of  this  study  demonstrate  the  ability  to  use  magnetic  fields  in  MERT  to  produce  highly  desirable 
dose  distributions. 


1.  Introduction 

The  current  trend  in  radiation  therapy  towards  conformal  techniques  requires  the  ability  to 
deposit  high  dose  to  a  given  target  region,  while  sparing  surrounding  normal  tissue.  Depth- 
dose  distributions  approaching  the  ideal  of  a  step  function  are  available  from  the  use  of  protons, 
heavy  ions  and  exotic  particles  such  as  tt  mesons.  However,  at  present,  these  treatment 
modalities  are  of  greatly  restricted  use  due  to  the  extreme  expense  and  technical  complexity 
of  the  accelerator  systems.  In  contrast,  an  electron  linear  accelerator  may  cost  a  factor  of  a 
hundred  less  than  a  proton  accelerator,  and  be  orders  of  magnitude  more  compact;  electron 
accelerators  are  thus  found  in  most  radiation  therapy  centres.  Furthermore,  while  photon  beam 
intensity-modulated  radiation  therapy  (IMRT)  allows  effectively  arbitrary  dose  distributions, 
the  requirement  of  several  centimetres  of  tissue  for  electron  build-up  makes  IMRT  of  shallow 
tumours  difficult.  Thus,  at  present,  electron  beams  hold  many  advantages  for  the  treatment  of 
shallow  tumours  such  as  those  in  the  breast  and  head  and  neck. 

Additionally,  a  recent  Monte  Carlo  investigation  (Ma  et  al  2000b)  into  modulated  electron 
beam  radiation  therapy  (MERT)  has  shown  that  the  use  of  non-uniform  intensity  profiles  and 
multiple  electron  beam  energies  will  allow  highly  conformal  dose  delivery  to  shallow  targets. 
Conformity  is  achieved  laterally  via  an  electron  multileaf  collimator  (Lee  et  al  2000a,  b) 
and  in  the  depth  direction  by  means  of  energy  and  intensity  modulation.  This  has  been 
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shown  to  be  particularly  useful  in  treatment  of  the  intact  breast.  MERT  would  avoid  both 
the  problems  of  dose  delivery  to  the  contralateral  breast  and  the  small,  very  high-dose  regions 
in  the  lung. 

Unfortunately,  while  electron  beams  exhibit  rapid  drops  in  depth-dose  distributions  relative 
to  photon  beams,  the  distance  between  the  end  of  the  therapeutically  useful  range  and  the 
practical  range  of  the  electrons  may  still  be  of  the  order  of  centimetres.  Thus,  while  MERT 
delivers  a  highly  conformal  dose  to  the  target,  because  of  the  low  density  of  the  lung  MERT 
exposes  much  of  the  lung  to  low  doses  of  radiation  (less  than  50%  of  the  overall  maximum.) 
The  biological  impact  of  this  low-dose  region  is  as  yet  unclear;  however,  it  is  likely  that 
the  clinical  outcome  will  improve  if  the  lung  dose  is  reduced.  It  is  thus  of  paramount 
importance  to  investigate  methods  that  may  allow  improved  depth-dose  distributions  from 
electron  beams. 

Numerous  studies  have  been  performed  to  examine  the  possibility  of  using  magnetic 
fields  to  tailor  the  dose  distributions  from  electron  beams.  Two  general  methods  for 
two  distinct  purposes  have  been  proposed:  magnetic  fields  directed  along  the  axis  of  the 
beam  (longitudinal)  for  lateral  confinement  of  penumbras  and  stability  in  the  presence  of 
inhomogeneities  (Weinhous  et  al  1985,  Bielajew  1993),  and  fields  directed  orthogonally  to 
the  beam  (transverse)  to  reduce  the  dose  penumbra  in  the  depth  direction.  In  this  work,  only 
the  effects  of  transverse  magnetic  fields  are  considered,  as  the  primary  goal  is  restriction  of 
electron  penetration.  Transverse  fields  were  suggested  by  Sempert  (1960)  and  early  Monte 
Carlo  studies  were  performed  by  Shih  (1975),  demonstrating  the  ability  to  create  a  sharp  depth- 
dose  drop-off.  Experimental  investigation  ensued,  including  work  by  Whitmire  et  al  (1977), 
who  conducted  experiments  with  22  and  28  MeV  beams  in  1  T  fields.  Nath  and  Schulz  (1978) 
performed  experimental  investigations  using  2.05  T  transverse  magnetic  fields  to  increase  the 
dose  peak  and  sharpen  the  dose  drop-off  of  an  electron  beam  of  less  than  50  MeV.  Further 
work  has  explored  the  effect  of  inhomogeneities  (Paliwal  et  al  1978,  Whitmire  et  al  1978). 
Recent  interest  in  Monte  Carlo  investigations  has  been  revived  by  the  work  of  Nardi  and 
Bamea  (1999),  using  the  ITS/ACCEPTM  system  to  demonstrate  the  use  of  modem  Monte 
Carlo  techniques  to  simulate  monoenergetic  electron  beams  in  a  magnetic  field  positioned  at 
a  fixed  depth. 

This  work  aims  to  use  Monte  Carlo  simulations  to  investigate  the  effects  of  transverse 
magnetic  fields  on  clinical  electron  beams,  and  to  examine  the  possibility  of  combining  this 
technique  with  MERT.  In  particular,  the  effect  of  different  field  strengths  and  field  locations  is 
examined.  The  changes  in  dose  distributions  are  quantified  with  a  focus  on  the  effect  on  surface 
dose  and  the  dose  fall-off.  A  modified  system  of  MERT  is  presented  in  which  magnetic  fields 
are  used  in  conjunction  with  intensity  modulation  to  generate  highly  conformal  dose  regions 
with  low  surface  dose  and  rapid  depth-dose  fall -off.  The  effect  of  heterogeneities  on  the  dose 
distribution  is  also  described. 


2.  Materials  and  methods 

2. 1.  Principles  and  theory 

When  an  electron  enters  a  magnetic  field,  it  immediately  begins  to  travel  in  a  curved  path  as 
dictated  by  the  Lorentz  force  and  eventually  exits  the  field  after  moving  along  a  semicircular 
path.  If  energy  loss  is  included,  then  the  particle  will  spiral  towards  smaller  and  smaller  radii. 
By  this  action,  an  electron  can  be  caused  to  deposit  much  of  its  energy  in  a  localized  region. 
In  particular,  the  electron  penetration  into  the  field  is  approximately  restricted  to  less  than  a 
distance  defined  by  the  radius  of  gyration,  r^. 
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In  the  limit  of  constant  energy  (small  step  size  A/),  in  the  presence  of  a  magnetic  field 
defined  by  5,  a  particle  direction  vector  Co  is  replaced  with  (Nelson  etal  1985) 

A/  ^-.A/ 

i;  =  5(i;o  *  +  [vo  —  B(i;o  *  j^)]  cos - vq  x  B  sin —  .  (1) 

As  Monte  Carlo  simulated  transport  occurs  in  short,  discrete  steps,  this  equation  is  particularly 
relevant  to  this  study. 

It  has  been  shown  by  Bielajew  (1993)  that  the  synchrotron  radiation  due  to  the  bending  of 
electron  paths  at  the  energies  of  interest  is  trivial:  of  the  order  of  0.1-1  eV  additional  stopping 
power.  Thus,  both  this  added  stopping  power  and  synchrotron  photon  dose  deposition  are 
taken  to  be  zero  in  this  study. 

One  may  generate  a  simple  mathematical  description  of  the  expected  effect  of  the  magnetic 
fields  on  depth-dose  distributions,  before  any  Monte  Carlo  simulations  are  performed.  Putting 
the  unperturbed  depth-dose  profile  in  a  vector  x  and  the  depth-dose  profile  in  a  magnetic  field 
in  a  vector  b,  there  exists  a  transformation  matrix  U  such  that 

Ujc  =  ^  (2) 

with  U  upper  triangular  with  columns  that  sum  to  unity.  The  effect  of  operating  on  x  with 
U  is  that  below  some  row  of  3c,  that  is,  beyond  a  given  depth,  there  is  a  ‘folding  back’  of 
energy  deposition,  such  that  some  portion  of  the  energy  that  would  otherwise  be  deposited 
at  greater  depths  will  now  be  deposited  earlier.  However,  because  the  overall  energy  of  the 
system  must  be  conserved,  the  columns  of  the  transformation  matrix  U  must  have  unit  sum. 
Additionally,  the  magnetic  field  cannot  cause  electrons  to  deposit  energy  further  downstream 
than  they  otherwise  would,  hence  the  upper-diagonal  nature  of  the  matrix.  It  will  be  seen 
the  effect  of  the  magnetic  field  is  a  function  of  the  incident  energy  (spectrum),  the  depth  the 
magnetic  field  begins  at,  and  the  strength  of  the  field.  In  principle,  if  such  a  function  could 
be  derived,  the  effects  of  magnetic  fields  could  be  computed  very  rapidly.  However,  for  the 
purposes  of  this  work,  such  analytical  formulations  are  considered  only  to  anticipate  the  form 
of  the  depth-dose  distributions  that  will  be  generated  by  Monte  Carlo  simulations. 

2.2.  Monte  Carlo  simulations 

Note  that  in  all  subsequent  sections,  for  brevity  and  clarity,  the  term  ‘beam’  will  be  used  to 
describe  the  incident  flux  of  electrons,  i.e.  the  electron  beam,  while  ‘field’  will  refer  to  the 
magnetic  field. 

Source-parameter  descriptions  of  electron  beams  were  obtained  according  to  a  procedure 
described  elsewhere  (Kapur  et  al  1998,  Jiang  et  al  2000,  Lee  et  al  2000c)  and  summarized 
here.  Electron  beam  simulations  of  a  Varian  CHnac  2100C  (Varian  Oncology  Systems,  Palo 
Alto,  CA)  were  performed  using  the  EGS4/BEAM  user  code,  (Nelson  et  al  1985,  Rogers 
et  al  1995)  using  vendor  supplied  geometries.  The  simulation  code  and  geometry  employed 
have  been  previously  shown  to  provide  agreement  with  measured  data  of  better  than  2%  in 
transverse  profiles  and  depth-dose  curves  for  the  three  nominal  energies  simulated:  6,  12  and 
20  MeV  (Kapur  et  al  1998).  Phase  space  files  were  obtained  above  the  last  scraper  of  electron 
applicators.  Source-parameter  files  were  then  obtained  by  using  the  BEAMDP  program  (Ma 
and  Rogers  1997)  and  used  in  all  future  dose  calculations.  The  four-source  model  used  in  this 
study  is  described  in  detail  in  Jiang  et  al  (2000).  The  source-parameter  files  have  also  been 
shown  to  agree  with  the  measured  data  to  better  than  2%  in  depth-dose  and  profiles  (Lee  et  al 
2000c). 

The  macro  package  for  electromagnetic  transport  developed  by  Bielajew  (1987,  1993) 
were  implemented  into  the  EGS4/MCDOSE  code  (Ma  et  al  2000a).  The  three-dimensional 
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Figure  1.  The  geometry  of  the  simulated  system.  The  magnetic  field  vector  is  orthogonal  to  the 
electron  beam  axis,  and  is  normal  to  the  plane  shown  in  the  isodose  curves.  As  shown,  the  direction 
of  electron  deflection  is  out  of  the  page.  The  magnetic  field  region  is  a  semi-infinite  slab,  i.e.  all 
points  below  a  given  depth  are  assigned  the  same  magnetic  field  vector. 


dose  distributions  were  examined  using  STATDOSE  (McGowan  et  al  1996).  Verification  of 
the  implementation  was  performed  by  computing  longitudinal  magnetic  fields  for  20  MeV 
electrons,  as  per  Bielajew  (1993)  and  15  MeV  monoenergetic  beam  in  fields  of  varying 
magnitude,  as  per  Nardi  and  Bamea  (1999)  using  ITS/ACCEPTM.  The  agreement  was  better 
than  2%  in  all  cases. 

Pencil  beams  were  simulated  by  restricting  the  particles  to  a  single  initial  track. 
Specifically,  a  monoenergetic  parallel  beam  was  used,  with  a  vanishing  aperture  size.  The 
voxel  geometry  employed  for  these  calculations  was  1x1x1  mm^.  The  narrow  beams 
were  simulated  by  using  the  source-parameter  file  for  a  6  x  6  cm^  applicator  and  transporting 
particles  through  a  3  x  3  cm^  Cerrobend  cut-out  before  entering  the  phantom.  The  10x10  cm^ 
open  fields  were  simulated  with  the  appropriate  source-parameter  file,  with  the  last  scraper 
included  in  the  MCDOSE  simulation.  The  voxels  for  the  narrow  and  broad  fields  were  2  mm 
in  the  depth  direction,  4  mm  in  the  axis  orthogonal  to  the  beam  and  magnetic  field  and  10  mm 
in  the  direction  of  the  magnetic  field.  When  the  effects  of  inhomogeneities  were  considered, 
specific  voxels  in  the  phantom  were  converted  to  either  air,  or  ICRP  standard  bone.  In  all  cases, 
the  nominal  source-to-surface  distance  (SSD),  defined  from  the  photon  target  to  the  phantom 
surface)  was  100  cm. 

In  all  cases,  the  magnetic  field  was  directed  orthogonally  to  the  electron  beam  axis,  and 
normal  to  the  plane  shown  in  the  isodose  curves.  The  geometry  is  shown  in  figure  1. 

All  Monte  Carlo  simulations  were  performed  on  various  elements  of  a  suite  of 
22  Pentium  Pro  (200  MHz)  and  10  Pentium  III  (450  MHz)  CPUs,  all  running  EGS4,  BEAM 
and  MCDOSE  under  the  Linux  operating  system.  All  simulations  utilized  an  electron  cut-off 
of  0.7  MeV  total  energy  and  a  photon  cut-off  of  10  keV,  below  which  all  remaining  energy 
was  deposited  on  the  spot.  PRESTA  extensions  were  employed  for  step  length  calculations 
(Bielajew  and  Rogers  1987).  In  general,  it  was  found  that  10  x  10^  electron  histories  were 
sufficient  to  provide  a  statistical  error  of  less  than  1%  (la).  The  overall  speed  of  calculation 
was  reduced  by  approximately  50%  relative  to  an  identical  simulation  with  no  magnetic  field 
in  place. 


2.5.  Optimization 

A  one-dimensional  optimization  routine  was  implemented  to  examine  the  feasibility  of 
combining  fields  such  that  the  narrow  dose  peaks  were  spread  out  to  several  centimetres. 
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Figure  2.  The  theoretically  derived  track  of  an  unscattered  electron  passing  through  water  in  the 
presence  of  a  1 .5  T  magnetic  field,  overlaid  on  isodose  curves  computed  using  the  Monte  Carlo 
method  for  a  pencil  beam:  (a)  6  MeV  with  magnetic  field  beginning  at  1  cm  depth,  (b)  12  MeV  at 
2  cm  depth,  (c)  20  MeV  at  3  cm  depth. 


The  dose  at  a  given  depth  in  the  combined  beam,  Act  was  taken  to  be  a  linear  combination  of 
the  dose  at  that  depth  of  the  individual  beams.  That  is,  for  n  different  field /beam  combinations 

n 

Aoi(z)  =  Xi  ■  Di{z).  (3) 

i=\ 

Each  beam  i  is  given  a  weight  x/,  which  is  determined  based  on  the  planned  dose  profile. 
To  determine  the  weights  x  =  (xi,  X2, . . . ,  x„),  depth-dose  curves  were  extracted  from  the 
10  X  10  cm^  fields  described  earlier  and  each  curve  placed  in  the  columns  of  a  matrix  A.  The 
desired  depth-dose  distribution  was  placed  in  a  column  vector  b  and  the  system 

Ax  ^  b  (4) 

was  solved  in  a  least-squares  sense  under  non-negativity  constraints.  The  method  of  solution 
utilized  is  described  in  detail  in  Lawson  and  Hanson  (1974). 

3.  Results  and  discussion 

The  basic  characteristics  of  electron  beams  in  magnetic  fields  are  explored  in  section  3.1 
and  section  3.2  by  examining  isodose  curves  and  depth-dose  curves.  This  demonstrates  the 


2952 


M  C  Lee  and  C-M  Ma 


Figure  3.  (a)  Isodose  curves  for  a  12  MeV  electron  beam  collimated  by  a  6  x  6  cm-  applicator  with 
a  3  X  3  cm"  Cerrobend  cut-out  in  the  presence  of  a  1 .5  T  magnetic  field.  Dotted  lines  represent  the 
central  axis  and  the  ray  passing  through  the  dose  maximum.  The  magnetic  field  is  directed  into  the 
page,  (b)  The  corresponding  depth-dose  curves  for  the  central  axis  and  the  ray  passing  through  the 
f^max  • 

general  principles  of  electron  beam  therapy  with  magnetic  fields,  and  provides  evidence  for  the 
viability  of  the  Monte  Carlo  transport  code.  Parameters  for  describing  the  depth-dose  curves 
are  presented  in  section  3.3  and  the  relevance  of  these  parameters  to  MERT  is  discussed.  An 
example  of  the  application  of  MERT  techniques  to  electron  beams  in  magnetic  fields  is  then 
given  in  section  3,4.  Practical  considerations  important  to  actual  application  of  magnetic  field 
MERT  are  presented  in  sections  3.5  and  3,6. 

5,7.  Pencil  beams  and  narrow  beams 

Monte  Carlo  calculations  were  performed  in  which  electrons  were  constrained  to  travel  down 
the  same  initial  axis,  i.e.  a  pencil  beam.  The  resultant  isodose  curves  are  shown  in  figure  2. 
These  figures  clearly  show  the  curvature  of  the  electron  beam  dose  distribution  generally  along 
the  expected  track.  By  superposition  of  many  pencil  beams  along  the  lateral  axis,  it  is  clear 
that  Bragg-peak-like  dose  peaks  can  be  constructed. 

Narrow  fields,  defined  by  3  x  3  cm^  Cerrobend  cut-outs  on  6  x  6  cm^  applicators,  were 
also  simulated.  What  is  apparent  is  that  there  is  a  significant  shift  of  the  dose  peak  off  the 
central  axis,  as  seen  in  figure  3.  These  data  are  presented  primarily  as  a  demonstration  of  the 
effects  of  the  magnetic  field,  i.e.  the  curvature  of  the  beam  and  the  formation  of  a  dose  peak. 
However,  it  is  difficult  to  make  any  definitive  statements  regarding  depth-dose  distributions 
because  of  the  lack  of  equilibrium,  as  discussed  below. 

5.2.  Broad  beams 

3.2  J.  Electronic  equilibrium.  In  the  3  x  3  cm^  narrow  beams  and  the  pencil  beams,  a  state 
of  lateral  disequilibrium  exists  such  that  along  any  given  axis,  many  times  more  particles  are 
removed  from  the  axis  by  the  magnetic  field  than  are  replaced,  as  seen  by  the  displacement  of 
the  dose  peak  in  figure  3,  In  the  limit  of  infinitely  broad  beams,  there  is  equilibrium  along  any 
ray  parallel  to  the  beam  axis.  This  is  also  true  in  the  central  region  of  broad,  finite  beams.  This 
can  be  observed  by  computing  the  integral  dose  along  a  central  axis  of  an  unperturbed  beam 
which  is  inherently  large  enough  for  equilibrium  to  exist.  If  lateral  equilibrium  exists,  then  the 
total  energy  deposited  along  this  axis  should  not  change  when  a  magnetic  field  is  applied. 
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Figure  4.  The  ratio  of  energy  deposited  in  the  central  axis  voxels  in  the  presence  of  a  magnetic 
field  to  the  energy  deposited  in  a  field-free  system.  This  ratio  is  shown  as  a  function  of  the  magnetic 
field  starting  depth:  (a)  1.5  T  magnetic  field,  (b)  3.0  T. 


Figure  5.  (a)  Lateral  dose  profile  for  a  20  MeV  beam  in  the  presence  of  a  1.5  T  magnetic  field 
beginning  at  3  cm  depth.  The  axis  being  profiled  is  orthogonal  to  both  the  beam  axis  and  the 
magnetic  field.  Shown  are  the  profiles  at  1.5,  2.5  and  3.5  cm  depths,  (b)  Depth-dose  curves  for  a 
12  MeV  beam  in  the  presence  of  a  1.5  T  magnetic  field  beginning  at  3  cm  depth.  Shown_^are  the 
central  axis,  and  axes  at  —4,  0, 2  and  3  cm,  where  the  positive  direction  is  defined  by  u  x  B  where 
V  is  the  original  beam  direction. 


The  Monte  Carlo  simulations  suggest  that  for  the  energies  investigated,  a  10  x  1 0  cm^  beam 
is  sufficiently  broad  to  overcome  the  disequilibrium  caused  by  the  magnetic  fields.  As  seen 
in  figure  4,  when  the  ratio  of  the  energy  deposited  along  the  central  axis  in  the  magnetic  field 
to  the  energy  deposited  in  the  zero-field  case  is  less  than  unity,  there  is  a  net  flow  of  electrons 
out  of  the  field.  This  is  in  fact  not  an  issue  of  lateral  equilibrium  in  the  standard  sense,  but  a 
danger  of  using  high-strength  magnetic  fields  at  shallow  depths.  What  is  occurring  is  that  the 
field  is  strong  enough  that  high-energy  electrons  are  directed  back  towards  the  surface,  and 
the  primary  beam  actually  exits  the  phantom  via  the  front  surface.  Clearly,  this  situation  is 
highly  undesirable  from  both  a  treatment  and  radiation  safety  standpoint,  and  it  is  important 
to  use  a  parameter  such  as  the  integral  dose  ratio  to  appreciate  the  minimum  safe  depth  for  a 
magnetic  field.  As  expected,  this  depth  is  shallower  at  higher  field  strengths,  as  more  electrons 
are  turned  back  towards  the  surface.  At  1.5  T,  the  minimum  safe  field  depths  are  0.4,  <  1  and 
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Figure  6.  The  80%  and  50%  isodose  curves  for  10  x  10  cm-  fields.  Both  the  zero-magnetic  field 
case  and  the  curves  with  a  1 .5  T  magnetic  field  directed  out  of  the  page  are  shown:  (a)  6  MeV  with 
magnetic  field  starting  at  1  cm  depth,  (b)  12  MeV  at  2  cm,  (c)  20  MeV  at  4  cm. 


1  cm  for  6,  12  and  20  MeV  respectively.  At  3.0  T,  the  depths  are  1,2  and  3  cm  at  6,  12  and 
20  MeV.  Note  that  some  points  have  values  slightly  greater  than  unity  (^1%  higher.)  This 
is  partially  due  to  the  photon  dose:  if  a  magnetic  field  causes  increased  electron  interactions 
near  the  surface,  the  photon  dose  curve  is  shifted  towards  the  surface.  Thus,  relatively  more 
photon  dose  will  be  deposited  in  the  phantom  in  this  case  than  in  the  field-free  case,  in  which 
the  photons  exiting  the  phantom  retain  more  energy.  This  effect  in  conjunction  with  statistical 
uncertainty  leads  to  the  slightly  higher  values. 

Once  the  problem  of  beam  loss  is  resolved,  the  original  issue  of  lateral  equilibrium  can  be 
considered.  Evidence  for  the  uniformity  of  the  field  is  provided  by  examining  the  lateral  dose 
profiles,  as  shown  in  figure  5(a).  These  profiles  are  taken  along  the  axis  orthogonal  to  both 
the  magnetic  field  and  original  beam  axis.  It  is  seen  that  within  the  7  cm  region  from  —4.5 
to  2,5  cm  (each  voxel  was  1  cm  in  width),  the  dose  was  homogeneous  to  within  5.8%  of  the 
average  dose  in  that  region.  Off-axis  depth-dose  curves  are  shown  in  figure  5(b).  Here,  the 
variations  in  the  depth-dose  curves  caused  by  moving  off  of  the  central  axis  are  comparable 
with  those  in  open  unperturbed  10  x  10  cm^  fields  while  within  the  equilibrium  region  (—4,5 
to  2.5  cm),  but  it  is  clear  that  there  is  a  rapid  drop  at  3  cm,  where  the  electron  dose  has  been 
shifted  strongly  in  the  positive  direction,  as  shown.  It  is  noted,  however,  that  even  at  that 
position,  there  still  exists  a  noticeable  dose  peak  and  dose  drop-off. 
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Figure  7.  Central  axis  depth-dose  curves  for  electron  beams  in  the  presence  of  a  1 .5  T  magnetic 
field.  For  each  energy,  the  depth-dose  curves  are  shown  for  different  magnetic  field  starting  depths. 
The  vertical  axis  represents  absolute  dose  in  arbitrary  units:  (a)  6  MeV,  (b)  12  MeV,  (c)  20  MeV. 


3.3,  Dose  profile  parameters 

Having  thus  established  the  electronic  equilibrium  of  these  systems,  it  is  possible  to  examine 
the  central-axis  dose  profiles.  Isodose  curves  for  10  x  10  cm^  fields  are  shown  in  figure  6. 
Profiles  along  the  central  axis  are  displayed  for  6,  12  and  20  MeV  electron  beams  with  the 
magnetic  field  positioned  at  several  depths  in  figure  7  (1.5  T)  and  figure  8  (3.0  T).  These  broad 
electron  fields  were  then  characterized  by  consideration  of  the  following  parameters,  defined 
graphically  in  figure  9: 

dg  =  the  depth  at  which  the  magnetic  field  begins 
^80/80  =  the  width  of  the  region  bounded  below  by  the  80%  dose  level 
^80/20  =  the  distance  between  the  deepest  80%  and  the  20%  dose  levels 
^p/eff  =  effective  practical  range,  found  by  extrapolation  of  the  drop-off  region 
SSF  =  surface  sparing  factor,  defined  by  the  ratio 

The  values  of  d^o/m^  d^o/20  and  SSF  for  selected  beam  and  field  combinations,  as  well 
as  the  no  magnetic  field  case,  are  shown  in  table  1  and  are  discussed  in  greater  detail  in  the 
following  sections. 

3.3.1.  Relationship  between  field  strength  and  dose  curves.  Monte  Carlo  simulations  were 
performed  in  which  the  magnetic  field  position  was  fixed  but  the  field  strength  was  varied. 
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Figure  8.  Central  axis  depth-dose  curves  for  electron  beams  in  the  presence  of  a  3.0  T  magnetic 
field.  For  each  enei^y,  the  depth-dose  curves  are  shown  for  different  magnetic  field  starting  depths. 
The  vertical  axis  represents  absolute  dose  in  arbitrary  units:  (a)  6  MeV,  (b)  12  MeV,  (c)  20  MeV. 


The  effect  of  the  changes  on  dso/20  and  SSF  are  shown  in  figure  10.  As  expected,  the 
distance  between  the  80%  and  20%  dose  falls  with  increasing  field  strength.  While  the  simple 
Lorentz  force  equation  may  suggest  a  linear  relationship,  due  to  energy  straggling  and  particle 
scattering,  a  more  complex  asymptotic  relationship  is  observed.  At  approximately  1.5  T,  the 
gains  from  increasing  the  magnetic  field  appear  to  be  minimal. 

The  SSF  shows  a  sigmoidal  behaviour,  with  a  period  of  rapid  increase  in  SSF  between 
0.75  and  1.5  T,  then  a  reduction  in  SSF  at  higher  field  strengths.  This  reduction  occurs  when 
the  radius  of  gyration  of  the  electron  is  too  small.  In  those  instances,  instead  of  travelling  in 
a  spiral,  the  electron  completes  an  approximate  half-circle  and  exits  the  field  at  some  angle, 
thus  depositing  dose  at  depths  shallower  than  the  magnetic  field. 

On  the  basis  of  these  results,  a  field  strength  of  1.5  T  was  chosen  as  the  primary  field  of 
interest.  It  is  noted  that  the  newest  generation  of  interventional  magnetic  resonance  imaging 
systems  produce  fields  of  this  magnitude.  Generally,  simulations  were  also  performed  at  3.0  T 
to  illustrate  in  detail  the  effects  of  changing  magnetic  field  strengths. 

3.3.2.  ^80/80 — therapeutic  range.  The  value  of  ^/go/so  was  taken  as  a  measure  of  the  portion 
of  the  beam  that  is  clinically  useful.  Values  for  selected  magnetic  field  positions  are  shown 
in  table  1.  This  region  is  extremely  narrow  in  all  cases,  with  a  maximum  value  of  1.7  cm  for 
a  20  MeV  beam  with  a  1.5  T  field  beginning  at  1  cm  depth.  As  a  result,  taken  as  is,  these 
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Figure  9.  The  parameters  used  to  describe  the  depth-dose  curves,  ^/go/so  is  the  distance  between 
the  two  80%  isodose  levels,  dso/20  is  the  distance  between  the  deepest  80%  and  the  20%  dose 
levels,  dB  is  the  depth  at  which  the  magnetic  field  begins  (it  is  taken  to  continue  to  infinite  depth), 
DS"  and  O,';-'  are  taken  as  approximate  measures  of  the  dose  maximum  and  minimum  along 
the  axis,  and  are  used  to  calculate  the  surface  sparing  factor  (SSF). 


Figure  10.  The  effects  of  changing  the  magnetic  field  strength  was  assessed  by  fixing  the  field 
position  and  incident  energy,  and  observing  the  effects  on  two  key  parameters:  (a)  the  value  of 
dm/2n  for  6  and  20  MeV  at  1.0  and  4.0  cm  depth  respectively;  (b)  the  ratio  of  peak  to  surface  dose, 
SSF,  for  6  and  20  MeV  at  1 .0  and  4.0  cm  depth  respectively. 


irradiation  geometries  are  of  limited  clinical  value.  However,  this  limitation  may  be  overcome 
using  the  modified  MERT  techniques  described  in  section  3.4. 


3.3.3.  dso/20 — dose  fall-off.  The  value  of  ^80/20.  the  width  of  the  80-20  ‘penumbra’  in  the 
depth  direction,  is  a  measure  of  the  ability  of  the  magnetic  field  to  improve  the  dose  drop-off, 
one  of  the  primary  goals  of  using  magnetic  fields.  It  was  found  that  with  a  3.0  T  magnetic  field, 
<^80/20  was  largely  invariant  as  the  depth  of  the  magnetic  field  was  changed,  as  seen  in  table  1. 
The  Monte  Carlo  results  show  an  average  ^^80/20  of  2.1  mm  at  6  MeV,  2.7  mm  at  12  MeV  and 
3.5  mm  at  20  MeV.  For  comparison,  note  that  without  a  magnetic  field,  the  d^o/io  values  are 
7.3  mm  for  6  MeV,  14  mm  for  12  MeV  and  27  mm  for  20  MeV.  This  indicates  a  reduction  in 
the  drop-off  distance  of  between  7 1%  for  6  MeV  to  87%  at  20  MeV. 
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Table  1.  Summary  of  beam  parameters  for  selected  electron  beams  at  1.5  and  3.0  T.  Parameters 
are  defined  in  section  3.3  and  are  given  along  the  central  axis. 


llB|i  =  1.5T  |lB[t  =  3.0T 


Energy 

(MeV) 

fi-field  starting 
depth  (cm) 

^/so/so 

(cm) 

^^8{)/:o 

(cm) 

SSF 

<^80/80 

(cm) 

^80/20 

(cm) 

SSF 

6 

0.0 

0.46t 

0,32^ 

l.50t 

0.24+ 

0.17^ 

1.28^ 

1.0 

0.49 

0.32 

2.35 

0.32 

0.25 

2.22 

12 

1.0 

0.84 

0.47 

2.31 

0.53t 

0.27t 

1.66^ 

2.0 

0.71 

0.51 

2.66 

0.46 

0.27 

2.37 

3.0 

0.67 

0.49 

2.25 

0.67 

0.28 

2.62 

4.0 

0.93 

0.38 

1.64 

0.73 

0.27 

2.01 

20 

1.0 

1.66^ 

0.66+ 

2.22+ 

0.91+ 

0.32+ 

1.61  + 

2.0 

1.39 

0.72 

2.55 

0.92+ 

0.35+ 

1.74+ 

3.0 

1.20 

0.74 

2.71 

0.82 

0.36 

2.12 

4.0 

1.06 

0.76 

2.57 

0.76 

0.40 

2.64 

5.0 

0.93 

0.72 

2.27 

0.62 

0.37 

2.63 

6.0 

0.96 

0.65 

1.88 

0.95 

0.34 

2.19 

7.0 

1.19 

0.58 

1.47 

0.84 

0.34 

1.81 

6 

None 

1.66 

0.73 

1.35 

1.66 

0.73 

1.35 

12 

None 

4.22 

1.40 

1.21 

4.22 

1.40 

1.21 

20 

None 

6.96 

2.70 

1.12 

6.96 

2.70 

1.12 

The  error  on  doses  (and  on  all  derived  values  shown  above)  was  in  all  cases  <1%  (lo)  of  the 
maximum  dose  in  the  volume. 

Lateral  disequilibrium.  Total  energy  deposited  along  central  axis  with  B-field  differs  from  the 
unperturbed  case  by  more  than  5%.  See  discussion  in  section  3.3. 


At  1 .5  T,  the  ^^80/20  for  the  6  MeV  beam  is  nearly  constant  at  3.3  mm.  However,  as  seen  in 
table  I,  there  is  a  noticeable  peak  in  the  ^/8o/20  values  for  12  and  20  MeV.  The  rise  in  ^^80/20  is 
due  to  the  increased  scattering  at  lower  electron  energies.  As  the  magnetic  fields  are  positioned 
deeper  in  the  phantom,  the  spectrum  of  the  electrons  caught  in  the  magnetic  field  is  shifted 
towards  lower  energies  and  the  increased  scattering  at  these  energies  results  in  a  scattering 
towards  greater  depths.  However,  as  the  magnetic  field  moves  beyond  this  peak  value  (2  cm  at 
12  MeV  and  4  cm  at  20  MeV),  the  energies  ('^8  MeV)  are  too  low  for  the  electrons  to  escape 
very  far  from  the  magnetically  defined  track. 

In  the  context  of  MERT,  regardless  of  how  many  fields  are  directed  along  the  same  port,  the 
d%o{20  of  the  combined  field  will  be  the  ^/8o/20  of  the  deepest  penetrating  field,  given  a  minimal 
contribution  from  the  bremsstrahlung  tail.  Thus,  before  any  combined  fields  are  simulated, 
there  is  an  a  priori  expectation  that  combined  fields  will  preserve  the  rapid  dose  fall-off  seen 
with  the  individual  fields. 


33.4,  Rp/eff— depth  of  maximum  penetration.  The  third  parameter  by  which  the  fields  were 
characterized  was  the  effective  practical  range.  This  is  similar  to  the  standard  dosimetric 
practical  range  Rp  in  that  it  is  defined  by  the  intersection  of  the  zero  dose  level  and  a  linear 
extrapolation  of  the  dose  fall-off  region.  Physically,  Rp/ea  represents  the  maximum  penetration 
of  electrons  in  the  medium  in  the  presence  of  the  magnetic  field.  Thus,  the  distance  Rp/ef^-de 
represents  the  distance  required  by  the  magnetic  field  to  turn  the  highest-energy  electrons  back 
towards  the  surface. 
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Figure  11.  The  surface  sparing  factor  (SSF),  defined  as  the  ratio  A  higher  SSF 

value  suggests  lower  skin  dose  for  a  given  target  dose.  The  SSF  is  shown  as  a  function  of  magnetic 
field  placement:  (a)  6  MeV,  (b)  12  MeV,  (c)  20  MeV. 


As  electron  energy  scales  approximately  linearly  with  depth  and  scales  linearly  with 
energy,  it  follows  that  the  distance  between  the  start  of  the  magnetic  field  and  the  effective 
practical  range,  Rp/eff-ds,  should  scale  linearly  with  dp.  Indeed,  this  was  the  case  as  is  shown 
in  table  1 . 


3,3,5.  SSF— surface  sparing  factor.  The  surface  sparing  factor  is  a  unitless  parameter  defined 
as  the  ratio  of  the  dose  maximum  to  the  dose  to  a  voxel  at  the  surface  (i.e.  the  average 
dose  delivered  to  the  first  2  mm  of  depth.)  However,  a  high  surface  sparing  factor  does 
not  immediately  suggest  that  a  high  dose  may  be  delivered  to  a  target  at  several  centimetres 
depth  while  avoiding  a  high  skin  dose.  As  described  above,  the  dose  peak  may  extend  only 
several  millimetres  in  the  depth  direction,  and  would  almost  certainly  not  be  clinically  useful 
as  a  single  field.  Combining  fields  inevitably  leads  to  a  reduced  SSF,  as  will  be  seen  later. 
Thus  it  is  crucial  that  the  beams  that  will  be  used  in  a  magnetic-MERT  system  exhibit  a  high 
starting  SSF,  before  combination  with  other  fields. 

As  seen  in  figure  1 1,  at  each  energy,  there  is  a  unique  starting  depth  for  the  magnetic  field 
that  results  in  a  maximal  SSF.  At  6  MeV,  the  SSF  increases  as  the  magnetic  field  is  moved 
deeper  until  a  maximum  occurs  between  0.8  and  0,9  cm  at  1.5  T  and  1.2  cm  at  3.0  T.  At  12  MeV, 
this  depth  was  1.0  cm  at  1.5  T,  or  2.0  cm  at  3.0  T,  and  at  20  MeV,  this  depth  was  3  cm  at  1.5  T 
and  4.0  cm  at  3.0  T.  The  existence  of  a  maximal  value  is  related  to  several  factors.  The  SSF 
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Figure  12.  Central-axis  depth-dose  curves  demonstrating  how  linear  combinations  of  10  x  10  cm^ 
electron  beams  in  magnetic  fields,  with  different  intensities  and  magnetic  field  placements,  can  be 
used  to  create  homogeneous  dose  distributions.  Note  that  both  the  SSF  ^^80/20 

are  significantly  improved  relative  to  a  field-free  situation:  (a)  1.5  T  only,  (b)  1.5  T  and  3.0  T, 
(c)  isodose  curves  for  the  1.5  T  optimization. 


drops  dramatically  when  the  magnetic  field  is  introduced  at  shallower  depths,  an  effect  related 
not  so  much  to  the  formation  of  a  dose  peak  as  the  deposition  of  energy  at  the  surface.  These 
cases  correspond  to  those  which,  as  discussed  earlier,  lack  electronic  equilibrium  as  electrons 
exit  the  beam  through  the  front  surface  of  the  phantom.  At  greater  depths,  electrons  are  at  a 
lower  energy  and  so  are  confined  to  smaller  regions,  thus  resulting  in  an  increased  dose  peak. 
This  is  offset  by  the  fact  that  at  deeper  depths  a  greater  percentage  of  the  energy  has  already 
been  deposited. 


3.4,  Combined  fields 

As  seen  when  examining  the  depth-dose  profiles,  the  depths  at  which  these  magnetically 
perturbed  electron  beams  are  therapeutically  useful  spans  a  limited  range.  Defining  a 
prescription  dose  of  80%  of  the  dose  maximum,  the  width  of  the  therapeutic  region  is  given  by 
z/go/80  and  can  be  less  than  1  cm,  as  seen  in  table  1  and  section  3.3,2.  Clinical  targets  are  often 
significantly  larger  and  so  a  method  is  required  to  create  wider  dose  peaks.  A  similar  problem 
is  encountered  in  proton  beam  therapy,  in  which  case  a  system  of  variable  beam  attenuation 
based  on  a  range  modulation  wheel  may  be  used  to  create  a  ‘spread-out  Bragg  peak’  of  arbitrary 
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Figure  13.  The  effect  of  an  air  cavity  heterogeneities  on  dose  distributions  with  and  without 
magnetic  fields.  Note  that  the  heterogeneities  do  not  change  the  shape  of  the  depth-dose  curve 
significantly  when  magnetic  fields  are  applied.  The  direction  of  electron  deflection  here  is  into 
the  page.  The  heterogeneity  extends  infinitely  into  and  out  of  the  page,  (a)  Isodose  curves 
without  magnetic  field,  (b)  with  1.5  T  magnetic  field  from  3  cm  (directly  after  the  air  cavity), 
(c)  corresponding  central-axis  depth-dose  curves.  The  vertical  axis  represents  absolute  dose  in 
arbitrary  units. 


width.  However,  for  electron  beams,  any  material  in  the  beam  path  would  significantly  distort 
the  field  shape  by  causing  additional  scattering  and  also  risk  the  production  of  bremsstrahlung 
photons  in  the  case  of  high-Z  materials.  However,  the  general  principle  of  beam  modulation 
remains  valid,  and  suggests  the  use  of  modulated  electron  radiation  therapy. 

Previous  discussions  of  MERT  have  considered  the  use  of  fields  with  non-uniform  energies 
and  intensity  distributions.  Here,  the  technique  is  extended  to  the  use  of  multiple  magnetic 
field  positions,  i.e.  a  single  port  may  include  a  20  MeV  field  with  a  magnetic  field  beginning 
at  3.0  cm  depth  for  a  certain  number  of  monitor  units,  and  then  another  20  MeV  field  with  the 
magnetic  field  located  at  a  depth  of  5.0  cm  delivered  down  the  same  axis,  etc.  It  is  instructive 
to  compare  this  methodology  with  the  more  familiar  method  used  in  proton  beam  therapy. 
In  that  system,  physical  blocks  of  different  thicknesses  are  used  to  shift  the  Bragg  peak  to 
different  depths,  and  intensity  modulation  is  provided  by  allowing  the  beam  to  dwell  on  a 
given  modulator  for  a  variable  time  period.  In  this  case,  the  only  difference  is  the  means 
of  moving  the  pseudo-Bragg  peak,  that  is,  by  moving  the  magnetic  field.  Beam  intensity  is 
modulated  simply  by  changing  the  number  of  monitor  units  delivered  with  the  magnetic  field  at 
a  given  position.  It  was  believed  that,  as  with  the  proton  beam,  useful  fields  can  be  constructed 
while  maintaining  the  desirable  depth-dose  drop-off  and  low  skin  dose. 
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Figure  14.  The  effect  of  a  bone  heterogeneity  on  dose  distributions  with  and  without  magnetic 
fields.  The  direction  of  electron  deflection  here  is  into  the  page.  The  heterogeneity  extends  infinitely 
into  and  out  of  the  page,  (a)  Isodose  curves  without  magnetic  field,  (b)  with  1 .5  T  magnetic  field 
from  3  cm,  (c)  corresponding  central-axis  depth-dose  curves.  The  vertical  axis  represents  relative 
dose. 


A  simple  one-dimensional  optimization  routine  was  utilized  to  generate  two  fields.  In 
one  case,  a  target  profile  was  generated  where  the  dose  would  rise  linearly  from  50%  at  the 
surface  (relative  to  the  maximum  along  the  profiled  axis)  to  100%  at  1.0  cm,  then  drop  to  zero 
at  6.2  cm.  The  optimized  solution  was  solved  numerically,  and  a  linear  combination  of  the 
1 .5  T  fields  was  generated  to  match  the  target  values  in  a  least-squares  sense.  The  results  are 
shown  in  figure  12(a).  Using  the  terminology  from  section  3.3,  the  width  of  the  treatment 
region,  was  found  to  be  5.0  cm,  followed  by  a  ^/80/20  drop-off  of  0.77  cm.  Between  the 
depths  of  1.6  cm  and  5,6  cm,  the  dose  was  homogeneous  to  within  ±5.4%  of  the  average  dose. 
The  corresponding  isodose  curves  for  this  situation  are  shown  in  figure  12(c). 

A  second  target  profile  was  also  constructed,  again  rising  linearly  from  50%  at  the  surface 
to  100%  at  1 .0  cm,  but  then  falling  to  zero  dose  at  a  depth  of  4.8  cm.  In  this  case,  both  1 .5  T  and 
3.0  T  magnetic  fields  were  permitted.  The  results  are  shown  in  figure  12(b),  with  a  Jgo/so  of 
3.7  cm  and  a  ^/80/20  of  0.43  cm.  Between  1.2  cm  and  4.6  cm  depth,  the  dose  was  homogeneous 
to  within  ±4.5%  of  the  average  dose. 

It  is  noted  that  to  achieve  the  same  target  coverage  as  in  figure  12(a)  without  magnetic 
fields,  a  20  MeV  beam  would  have  to  be  used.  Recall  that  the  ^/8o/20  for  a  single  20  MeV  beam 
with  a  magnetic  field  is  2,7  cm,  or  250%  of  the  ^80/20  for  the  1.5  T  field  case.  Additionally, 
the  surface  dose  would  be  higher  by  approximately  10%  of  the  dose  maximum.  Similarly, 
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for  figure  12(b),  a  12  MeV  beam  would  be  used,  with  a  d^o/20  of  nearly  10  times  that  of  the 
magnetic  field  case,  with  a  significantly  higher  surface  dose  as  well.  MERT  techniques  may 
improve  the  dose  homogeneity  within  the  target  region,  but  cannot  improve  the  depth-dose 
fall-off. 

The  small  dose  inhomogeneities,  in  particular  small  peaks  near  the  end  of  the  treatment 
field,  are  caused  by  the  limited  dataset.  The  optimization  was  restricted  to  setting  the  magnetic 
field  at  the  fixed  depths  used  in  this  study.  A  full-scale  optimization  routine  would  allow  the 
standard  non-uniform  intensities,  in  addition  to  computing  optimal  magnetic  field  placements 
for  each  segment,  with  the  ability  to  choose  the  field  placement  along  a  continuum  or  among 
a  closely  spaced  set  of  discrete  depths. 

Certainly,  a  real  magnetic  field/beam  profile  optimization  would  have  to  occur  at  a 
minimum  in  thus  accounting  for  the  lateral  shift  seen  in  figure  5(a).  In  an  engineering 
sense,  reversing  the  magnetic  field  direction  is  likely  to  be  trivial,  and  this  possibility  must  also 
be  accounted  for.  The  mathematical  and  physical  nature  of  this  optimization  scheme  is  beyond 
the  scope  of  this  work,  but  further  discussion  of  practical  issues  involved  in  optimization  and 
treatment  planning  is  presented  in  section  3.6.2. 


3.5.  Inhomogeneities 

The  presence  of  shallow  heterogeneities  located  between  the  surface  and  the  magnetic  field 
was  also  investigated.  Slabs  of  bone  and  air  cavities  were  simulated  by  replacing  the  water 
of  some  of  the  voxels  with  air  or  ICRP  standard  bone.  The  goal  was  to  observe  the  effect  of 
the  changes  in  electron  energy  spectra  caused  by  passage  through  heterogeneous  media  on  the 
parameters  discussed  in  section  3.3. 

The  isodose  curves  with  an  air  cavity,  both  with  and  without  magnetic  fields,  are  shown 
in  figures  13(a)  and  13(b).  The  corresponding  depth-dose  curves  are  shown  in  figure  13(c). 
The  results  of  the  simulations  with  a  3  cm  slab  of  bone  at  the  surface  are  shown  in  figure  14. 

Note  that  the  goal  of  these  simulations  differs  from  previous  experimental  work,  such  as 
that  of  Paliwal  et  al  (1978).  In  that  work,  the  authors  used  a  magnetic  field  gradient  such 
that  the  dose  distribution  with  the  inhomogeneity  did  not  differ  significantly  from  the  dose 
distribution  without  the  inhomogeneity.  Here,  the  dose  distributions  differ  radically  between 
the  heterogeneous  media  with  the  magnetic  field  and  the  homogeneous  media  with  no  field. 
Instead,  the  primary  point  of  interest  is  that  the  general /orm  of  the  depth-dose  distributions  is 
preserved;  that  is,  a  dose  peak  followed  by  a  sharp  dose  drop-off  with  a  favourable  SSE  Given 
this,  then  even  in  the  presence  of  inhomogeneities  treatment  fields  can  be  constructed  by  using 
the  modulation  techniques  described  earlier. 

As  shown  in  figure  13(c),  with  an  air  cavity  present,  there  is  a  trivial  ( 1  mm)  decrease  in  the 
dmjio  and  a  decrease  in  the  SSF  from  2.7  to  2.2.  With  a  bone  slab,  there  was  no  change  in  ^80/20 
and  a  drop  in  SSF  from  2,7  to  2.5.  Thus,  even  if  the  beam  passes  through  inhomogeneities  before 
entering  the  magnetic  field,  the  advantages  of  magnetic  fields  are  still  preserved.  Certainly, 
with  an  inhomogeneity  present,  additional  fields  may  be  necessary  to  cover  the  entire  target 
region.  However,  it  has  been  shown  elsewhere  (Ma  et  al  2000b)  that  the  non-uniform  intensity 
profiles  used  in  MERT  can  provide  highly  conformal  coverage  of  a  target  in  the  presence 
of  heterogeneities.  Based  on  these  results,  it  is  concluded  that  MERT  combined  with  the 
techniques  shown  in  the  previous  section,  will  exhibit  the  conformal  coverage  of  jviert  with 
the  added  ability  to  reduce  surface  dose  and  dose  to  distal  structures,  even  in  inhomogeneous 
media. 
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Figure  15.  The  effect  of  finite  magnetic  field  gradients  on  central-axis  depth-dose  curves.  Gradients 
are  displayed  as  magnetic  field  strength  as  a  function  of  depth:  (a)  gradient  1 ,  (b)  gradient  2.  (c)  The 
corresponding  depth-dose  curves  with  the  infinite  gradient  included  for  comparison.  The  vertical 
axis  represents  relative  dose. 


3.6.  Practical  considerations 

3.6.1.  Field  gradients.  For  the  simulations  described  thus  far,  the  structure  of  the  magnetic 
fields  have  been  taken  to  be  a  step  function,  i.e.  the  field  gradient  is  infinite.  However,  this 
cannot  be  the  case  in  any  real  system.  Similarly,  the  field  strength  outside  the  primary  region 
is  typically  non-zero,  as  low-strength  fringe  fields  will  inevitably  fill  the  remaining  space. 
Nardi  and  Bamea  (1999)  demonstrated  that  for  a  15  MeV  monoenergetic  beam  in  a  3.0  T 
field,  the  difference  between  the  infinite  and  finite  gradient  (d||B||/dz  ^  1.25  T  cm“*)  was 
not  significant.  The  difference  in  Jgo/so  was  1  mm  while  the  SSF  dropped  from  2,43  to  2.27. 
As  demonstrated  above,  ^/so/so  can  be  made  largely  irrelevant  by  combining  fields.  Thus,  this 
section  focuses  on  the  effects  of  finite  gradients  on  ^30/20  and  SSF. 

To  explore  the  effects  of  finite  gradients  and  fringe  fields,  two  magnetic  fields  were 
constructed  with  gradients  shown  in  figures  15(a)  and  15(b).  Considering  the  semil inear 
region  centred  around  4  cm,  gradient  1  rises  with  d||B||/dz  0.4  T  cm”^  while  gradient  2 

rises  at  ~0.2  T  cm"^  The  resulting  central  axis  dose  profiles  are  shown  in  figure  15(c), 
compared  with  the  20  MeV  4  cm  field  case  (4  cm  corresponds  to  the  50%  field  level  in  the 
gradients.)  It  can  be  seen  that  the  ^30/20  drop-off  increases  by  less  than  0.5  mm  with  gradient  1 
(the  steeper  gradient),  and  less  than  1.5  mm  with  gradient  2.  While  this  is  considered  to 
be  an  adverse  effect,  these  ^^30/20  values  still  represent  70.3%  and  67.0%  improvements  over 
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the  zero-field  case.  The  peak  dose  value  was  also  seen  to  be  diminished  while  surface  dose 
remained  relatively  constant.  The  SSF  when  the  field  was  defined  by  a  step  function  was  2.57, 
With  gradient  1  this  was  reduced  to  2.31,  and  with  gradient  2  it  was  reduced  to  2.26.  Again, 
this  is  undesirable,  but  the  SSF  is  still  more  than  twice  that  of  the  20  MeV  zero-field  system. 

Thus,  it  is  concluded  that  while  sharper  gradients  are  clearly  superior  to  slower  gradients, 
even  an  unrealistically  ‘poor’  field  of  only  0.2  T  cm”^  still  presents  significant  advantages 
over  the  zero-field  case. 


3,6.2.  Other  considerations.  It  is  clear  that  the  largest  hurdles  in  the  implementation  of 
magnetic  fields  are  not  theoretical  issues,  but  rather  are  of  a  much  more  practical  nature. 
In  particular,  the  effect  of  fringe  fields  on  the  electronic  equipment  of  the  accelerator  and  the 
stability  of  the  electron  beam  as  it  travels  through  the  treatment  head.  This  must  be  investigated 
thoroughly  before  any  clinical  implementation  can  occur. 

Additionally,  the  use  of  magnetic  fields  at  multiple  depths  for  different  segments  of  an 
intensity  modulated  field  requires  high-precision  computer  control.  It  is  feasible  that  linking 
the  electron-specific  MLC  to  be  used  in  MERT  (Lee  etal  2000a,  b)  and  a  magnetic  field  system 
to  the  same  controller  will  allow  such  delivery  to  occur.  Thus,  for  example,  a  single  ‘step’  of  a 
step-and-shoot  type  algorithm  may  represent  the  motion  of  the  leaves  to  a  new  pattern  as  well 
as  motion  of  the  magnetic  field  to  a  new  position.  This  would  of  course  necessitate  significant 
changes  in  the  entire  treatment  planning  system.  For  example,  a  beamlet-based  beam  profile 
optimization  routine  must  have  precalculated  beamlets  not  only  for  different  energies,  but 
different  field  positions  (see  section  3.4).  The  leaf  sequence  algorithm  must  be  extended  to 
also  compute  the  necessary  movement  of  the  magnetic  field.  As  has  been  shown  elsewhere 
(Lee  et  al  2000c),  the  background  dose  can  rise  significantly  if  too  many  field  segments  are 
used,  so  the  treatment  planning  system  must  balance  the  background  dose  with  trying  to  achieve 
maximum  conformity.  Finally,  the  issue  of  magnetic  field  verification  is  created,  which  does 
not  lend  itself  to  as  straightforward  a  solution  as  leaf  sequence  verification. 

A  simplification  of  the  system  could  conceivably  utilize  magnetic  fields  simply  as  a  way 
of  restricting  dose  to  an  organ  at  risk.  For  example,  in  an  anterior-posterior  directed  MERT- 
based  breast  treatment  (Ma  et  al  2000b),  the  magnetic  field  could  be  aligned  along  the  base  of 
the  chest  wall.  This  will  not  improve  homogeneity  in  the  breast  or  reduce  skin  dose,  but  will 
restrict  dose  to  the  lung.  The  advantages  of  this  would  have  to  be  compared  with  the  potential 
dose  boost  to  the  distal  end  of  the  target  structures. 


4.  Conclusion 

It  has  been  demonstrated  that  transverse  magnetic  fields  can  be  used  effectively  to  create  dose 
distributions  with  sharp  depth-dose  fall-offs  and  a  favourable  ratio  of  peak-to- surface  dose.  The 
depth-dose  distribution  from  such  a  field  was  qualitatively  observed  to  resemble  the  Bragg  peak 
seen  in  proton  therapy.  The  position  and  size  of  this  peak  varied  as  function  of  electron  energy, 
field  strength  and  field  position.  Generally,  this  peak  exhibited  only  a  very  narrow  therapeutic 
range,  and  so  a  method  of  broadening  this  peak  was  necessary.  It  was  demonstrated  that  a 
modified  form  of  MERT  could  be  employed  to  generate  a  clinically  useful  field.  The  energy 
and  intensity  of  the  electron  beam,  as  well  as  the  position  of  the  magnetic  field  were  varied, 
thus  creating  a  homogeneous  dose  distribution  several  centimetres  in  depth,  while  retaining 
the  desired  dose  fall-off  and  low  surface  dose.  This  could  be  accomplished  with  constant 
magnetic  field  strength  or  variable  field  strength.  It  was  further  shown  that  the  primary  effect 
of  heterogeneities  is  to  shift  the  position  of  the  dose  peak,  but  to  have  minimal  effect  on  the  dose 
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fall-off  distance.  The  effect  of  magnetic  field  gradients  was  seen  to  adversely  affect  the  drop-off 
and  surface  sparing,  but  only  moderately.  Nonetheless,  it  is  clear  that  the  sharper  gradients  are 
preferred.  With  recent  improvements  in  the  development  of  compact  high-strength  magnetic 
fields,  higher- Tc  superconductors,  the  development  of  variance  reduction  techniques  and  rapid 
Monte  Carlo  electron  dose  calculation  systems  and  faster  computer  systems,  it  is  anticipated 
that  the  use  of  magnetic  fields  will  rapidly  become  practically  feasible.  Thus  it  is  concluded 
that  when  combined  with  the  emerging  modality  of  MERT,  magnetic  fields  allow  significant 
improvements  in  dose  conformity  in  the  depth  direction. 
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